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Abstract 


-DfiEA  has,  at  present,  two  libraries  containing  subroutines  for  calculating  splines: 
IMSL  and  BSPLIN.  A  new  library  has  been  developed  to  supplement  the  IMSL  and  BSPLIN 
routines  in  the  realm  of  smoothing  splines.  It  is  not  self-contained,  making  frequent  use 
of  subroutines  from  the  BSPLIN  library. 

The  new  subroutines  offer  several  advantages  over  the  smoothing  spline 
subroutines  in  the  IMSL  and  BSPLIN  libraries’, 

1 )  The  order  of  the  spline  may  be  picked  by  the  userj 

2)  The  second  derivative  of  the  spline  is  not  constrained  to  be  zero  at  its 
end-points^’ 

3)  The  user  of  the  new  subroutines  has  freedom  to  choose  the  number  and 
positions  of  the  knots  of  the  splinej  0,r  ■L 

4)  The  new  subroutines  have,  as  input,  an  extra  set  of  weights,  6,,  i*1,N, 
which  control  the  stiffness  of  the  spline  between  each  pair  of  knots. 

The  new  subroutines  were  initially  developed  for  use  in  ship  hull  approximation  for 
the  calculation  of  boundary  layer  growth  on  the  hull.  For  this  calculation  one  needs 
splines  whose  second  derivatives  are  very  well  behaved.  The  additional  control  afforded 
by  the  new  subroutines  makes  them  far  more  suitable  for  this  application  than  any  of  the 
subroutines  currently  available  in  either  the  IMSL  or  BSPLIN  libraries. 


Resume 


L'ERDA  possede  maintenant  deux  bibliotheques  contentant  des 
sous-programmes  pour  calculer  des  splines  :  IMSL  et  BSPLIN.  Une  nouvelle 
bibliotheque  a  ete  mise  sur  pied'  pour  completer  les  programmes  ISML  et 
BSPLIN  dans  le  domaine  des  splines  de  lissage.  Elle  n'est  pas  autonome, 
faisant  souvent  appel  a  des  sous-programmes  de  la  bibliotheque  BSPLIN. 

Les  nouveaux  sous-programmes  off rent  plusieurs  avantages  par 
rapport  aux  sous-programmes  de  splines  de  lissage  des  bibliotheques  IMSL 
et  BSPLIN. 


(1)  L'utilisateur  peut  choisir  le  degrA  de  la  spline. 

(2)  La  deuxieme  derivee  de  la  spline  n'est  pas  forcement  nulle 
a  ses  points  extremes. 

(3)  L’utilisateur  des  nouveaux  sous-programmes  peut  choisir  le 
nombre  et  le  lieu  des  noeuds  de  la  spline. 

(4)  Les  nouveaux  sous-programmes  acceptent  en  entree  un 
ensemble  supplAmentaire  de  coefficients  de  pond&ration  4^ 

i*l,  N,  qui  determinent  la  raideur  de  la  spline  entre  deux 
noeuds . 


Les  nouveaux  sous-programmes  ont  intialement  ete  mis  au  point 
pour  1* approximation  des  coques  de  navire,  notamment  pour  le  calcul  de  la 
croissance  des  couches  limites  sur  les  coques.  Pour  ce  dernier  calcul, 
il  faut  utiliser  des  splines  dont  la  deuxifcme  dArivAe  est  parfaitement 
dAfinie.  Par  le  contrdle  accru  qu'ils  offrent,  les  nouveaux 
sous-programmes  conviennent  beaucoup  mieux  k  cette  application  que  tous 
les  sous-prograimnes  existants  des  bibliotheques  IMSL  ou  BSPLIN. 
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NOTATION 


Bnk  "  The  ntt1  B-sp line  of  order  k  (Section  2.1) 

Dn^  "  The  m,h  divided  difference  operator  (Section  3) 

en  -  Error  of  the  n”1  data  point  (Section  2.1 ) 

f(x)  -  The  spline  function  (Section  2.1) 

F  -  The  smoothing  functional  as  used  by  Reinsch  and  de  Boor  (Section  2.1) 

F*  -  The  smoothing  functional  as  used  in  BSMTH  (Section  2.1) 

-  See  equation  (3.9) 

G  -  pX2  ♦  (1  -p)F  (Section  2.1) 

G*  -  pX2  ♦  (1-p)F*  (Section  2.1) 

k  -  The  order  of  the  spline  (Section  2.1 ) 

m  -  The  derivative  of  the  spline  function  used  as  a  smoothing  criterion 

(Section  2.1 ) 

N  -  The  number  of  B-splines  used  in  the  spline  (Section  2.1) 

Nk  -  The  number  of  knots  (Section  2.1) 

Np  -  The  number  of  data  points  (Section  2.1) 

p  -  Parameter  which  balances  the  relative  values  of  F  and  X2  (Section  2.1) 

phi  -  Value  of  p  used  in  the  Iteration  for  p  In  BSM.TH  (Section  2.5) 

p,0  -  Value  of  p  used  In  the  iteration  for  p  in  BSMTH  (Section  2.5) 

phl  -  The  value  of  p  after  the  n*  iteration  for  p  in  BSMTH  (Section  2.5) 

S  *  The  value  of  the  X2  input  by  the  user  (Section  2.5) 

sn  *  The  arc  length  to  the  0th  data  point  (Section  5) 

vi 

>  - 


V, 


-  See  equation  (2.10) 


v*„  -  See  equation  (2.19) 

wn  -  Error  weights  defined  in  equation  (3.2) 

xn  -  Abscissa  of  the  nth  data  point  (Section  2.1) 

Y2  -  See  equation  (2.1 1 ) 

Y*2  -  See  equation  (2.20) 

yn  -  Ordinate  of  the  nth  data  point  (Section  2.1) 

y,'n  -  See  equation  (2.21) 

a  -  Parameter  which  balances  the  relative  values  of  F  and  X2  (Section  2.1) 

|9n  -  The  nth  spline  coefficient  (Section  2.1) 

(9*n  -  Approximation  to  /3n  used  for  numerically  stable  determination  of  the  X2  of  a 

spline  (Section  2.4) 

6n  -  The  stiffness  weight  corresponding  to  the  interval  between  the  (n+k-1)tt1  and 

the  (n+k)th  knot  (Section  2.1) 

6nj  -  The  Kronecker  delta  (Section  3) 

£n  -  The  actual  error  of  the  n,h  data  point  (Section  3) 

O'  -  See  equation  (3.3) 

X2  -  The  chi-square  of  the  spline  (Section  1) 

X2h,  -  The  chi-square  of  the  spline  corresponding  to  the  p-value  phi  (Section  2.5) 

X2,0  -  The  chi-square  of  the  spline  corresponding  to  the  p-value  p,0  (Section  2.5) 

X2,  -  The  chi-square  of  the  spline  corresponding  to  the  p-value  pn  (Section  2.5) 
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1  INTRODUCTION 


DREA  has,  at  present,  two  libraries  containing  subroutines  for  calculating  splines: 
BSPLIN:  and  IMSL2.  The  library  of  subroutines  presented  here  is  intended  to  supplement 
the  previous  two.  It  is  not  self-contained,  making  frequent  use  of  BSPUN  subroutines. 

The  subroutines  presented  here  were  developed  because  the  smoothing  spline 
routines  available  in  IMSL  and  BSPLIN  were  found  inadequate  for  smoothing  data  digitized 
from  offset  diagrams  of  ship  hulls.  The  spline  representations  of  the  hulls  were  to  be 
used  in  the  calculation  of  hull  boundary  layer  growth.  For  this  application,  it  is  necessary 
to  have  a  spline  representation  of  the  hull  whose  second  derivatives  are  very  well 
behaved.  The  second  derivatives  of  the  hull  representation  cause  accelerations  in  the 
fluid  flow  around  the  hull  which  in  turn  cause  changes  in  the  boundary  layer  growth.  It 
was  found  that  the  spline  subroutines  in  the  IMSL  and  BSPLIN  libraries  could  not  be 
controlled  sufficiently  well  that  the  boundary  layer  calculations  would  be  unaffected  by 
splining  errors.  In  particular,  the  splines  were  unable  to  turn  sharp  corners  (near  the 
bilge,  for  example)  sufficiently  rapidly  without  either  cutting  the  corner  or  having 
'wiggles'  on  each  side  of  the  corner.  Either  result  induced  large  errors  in  the  second 
derivatives  of  the  spline,  the  former  underestimating  the  magnitudes  of  the  second 
derivatives,  the  latter  overestimating  them.  It  was  therefore  necessary  to  develop  new 
subroutines  providing  greater  control  over  the  splines  and  their  derivatives. 

The  most  fundamental  subroutine  in  the  new  library  is  BSMTH.  It  is  very  similar  in 
function  to  the  IMSL  subroutine  ICSSCU  (this  is  an  implementation  of  a  program  originally 
written  by  Reinsch3)  and  the  BSPLIN  subroutine  SMOOTH:  given  the  X2  of  the  spline  curve 
with  respect  to  given  data,  a  smooth  spline  approximating  the  data  is  determined  by 
minimizing  a  functional  which  measures  the  'lack  of  smoothness'  of  the  spline.  BSMTH, 
however,  offers  several  advantages  over  the  other  two  subroutines. 

1 )  The  order  of  the  spline  may  be  picked  by  the  user.  SMOOTH  and  ICSSCU 
are  cubic  splines  only. 

2)  SMOOTH  and  ICSSCU  constrain  the  second  derivative  of  the  spline  to  be 
zero  at  its  end-points.  BSMTH  imposes  no  such  constraint. 

3)  The  user  of  BSMTH  has  freedom  to  choose  the  number  and  positions  of  the 
knots  of  the  spline.  SMOOTH  and  ICSSCU  require  exactly  one  knot  at  each 
data  point.  The  freedom  to  choose  the  knots  allows  much  greater  control 
of  the  spline. 

When  splining  in  two  dimensions,  control  of  the  knots  has  additional 
consequences.  For  efficient  approximation  of  two-dimensional  data,  the 
knots  must  form  a  rectangular  lattice  (see  Reference  1,  chapter  17,  for 
example).  ICSSCU  and  SMOOTH  then  require  the  data  points  to  be  in  a 
rectangular  lattice.  With  BSMTH  this  is  no  longer  necessary. 

4)  BSMTH  has,  as  input,  an  extra  set  of  weights,  S|t  i=1,N,  which  control  the 
stiffness  of  the  spline  between  each  pair  of  knots.  If  the  spline  is 
required  to  be  flat  in  some  region,  then  the  appropriate  5,  is  increased.  If 
the  spline  is  to  bend  sharply  in  a  different  region,  the  appropriate  6,  is 
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=  1,  if  P  is  to  be  recalculated. 

Via  COMMON  /  PLIMS  / 

PMIN  =  Minimum  allowed  value  of  p  (See  Section  2.5).  Default  is  1.0E-03 

PMAX  =  Maximum  allowed  value  of  p.  Default  is  1.0E+03. 

Via  COMMON  /NODFLT/ 

IMAX  :  2*  I  MAX  is  the  maximum  number  of  divided  differences  allowed  to 
find  the  error  in  function  PRERR  (See  Section  A. 4).  Default  value 
is  5. 

SMFACT :  The  value  of  the  smoothing  parameter  used  by  BSMTH  may  be 
adjusted  by  using  a  value  of  SMFACT  not  equal  to  1.  The  smoothing 
parameter  used  is,  S  =  SMFACT*NPT*PRERR**2.  The  default  value 


Via  COMMON  /INTEXP/ 

JDER  :  The  value  of  JDER  used  by  BSMTH.  The  integral  of  the  square  of 
the  JDERth  derivative  of  the  spline  is  minimized  (subject  to  the 
constraint  that  XSQ  =  S).  If  smooth  curves  are  desired  a  value  of 
JDER  =  2  is  appropriate.  JDER  should  be  non-negative  and  less  than 
K.  The  default  value  is  2. 

DEFAULTS 

If  IER  s  0  on  input  then 
JDER  =  2 
SMFACT  *  1 .0 
IMAX  =  5 

T(l)  =  (l-K)/(N-K+1 ),  1=1, NKT  i.e.  knots  are  uniformly  distributed  in  (0,1) 

If  IER  =  1  on  input,  then  the  vales  for  JDER,  SMFACT,  IMAX  and  T(l)  must  be 
input  by  the  user  via  the  COMMON  blocks  /NODFLT/  and  /INTEXP/. 

OUTPUT 

IER  =  0,  Calculation  has  been  successful 
=  1 ,  If  JDER  >  K  -  1 

=  2,  If  NKT1  1  N  +  K  *  max(0,K-2*JDER) 

=  3,  If  IWK.  <  max(NKT1,K**2) 

=  4,  If  more  than  30  iterations  are  required  to  find  the  correct  value 
for  p  in  BSMTH  when  splining  the  data  point  abscissae.  Indicates 
numerical  difficulties  in  the  solution  of  the  linear  system 
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Section  6 


Appendix  A 
USER'S  GUIDES 


Concise  guides  for  the  use  of  the  spline  subroutines  are  now  given.  The 
subroutines  are  listed  alphabetically. 


A.1  BSMCRV  :  User's  Guide 

SUBROUTINE  BSMCRV(NPT,X,Y,E,N,K.NKT1  ,T,WTI,BC0EFX,BC0EFY,R,IWK.,WK.,ARCL,3,IER) 

PURPOSE:  Given  date  points  (X(I),Y(I)),  1=1, NPT  BSMCRV  finds  a  smooth  curve 

approximating  them  by  spiining  the  abscissae  and  ordinates  separately  with  respect 
to  the  arc-length  along  the  spline.  The  arc  length  at  each  point  is  approximated 
from  the  distances  between  the  points.  BSMTH  is  used  to  spline  the  abscissae  and 
the  ordinates.  The  function  PRERR  is  used  to  determine  the  smoothing  parameter 
and  the  subroutine  WTIBEG  is  used  to  determine  the  stiffness  weights. 

LANGUAGE: FORTRAN 

USAGE:  EXECUTE  mainpgm,BSPLIN:HLLYSP/UB,BSPLIN:BSPLIN/LIB 

CALLS  subroutines  BSMTH,  PRERR,  WTIBEG 

INPUT 


NPT 

X 

Y 

E 

N 

K 

NKT1 


:  The  number  of  data  points. 

•  An  array  of  length  NPT  containing  the  data  point  abscissae  in 
ascending  order. 

:  An  array  of  length  NPT  containing  the  data  point  ordinates. 

.  The  errors  of  the  data  points.  The  smaller  the  error  the  closer  the 
spline  will  come  to  that  point. 

:  The  number  of  B-splines  used  to  represent  the  spline. 

:  The  order  of  the  spline. 

=  N  +  K  +  max(Q,K-2*JDER)  (see  below  for  a  definition  of  JDER) 


T  :  An  array  of  length  NKT1  the  first  N+K  elements  of  which  contain  the 

knot  sequence  (in  ascending  order).  The  remaining  array  elements 
are  used  in  subroutine  SETUPP. 

IER  =  0,  If  JDER,  T,  WTI  and  the  first  N*K  elements  of  R  are  as  on  the 
previous  call  to  BSMTH  (this  means  that  the  matrix  P  need  not  be 
recalculated). 
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5  A  PARAMETRIC  SMOOTHING  SPLINE 

It  is  often  desired  to  approximate  data  by  a  smooth  curve  which  is  not  necessarily 
a  function.  The  spline  approximation  must  be  parametrized  in  some  way.  The  choice  of 
the  parametriztion  is  important  (see  Reference  1,  pp.316).  It  has  been  shown  that  any 
approximation  of  the  arc  length  of  the  curve  provides  a  good  parametrization.  It  is 
usually  sufficient  to  approximate  the  arc  length  from  the  distance  between  data  points. 
The  parametric  spline  is  then  calculated  as  follows. 

1 )  Calculate  the  parameter  s,  at  each  data  point  by 

sn  =  sn.,  +  (Cyn  -  yn--,)2  +  (y*  -  yn-i>2>H  t5-1 ) 

2)  Spline  each  of  the  data  sets  {(Sn.xJ.nsl.N)  and  {(sn,yn),n=1,N}. 

This  is  perfomed  in  the  subroutine  BSMCRV,  which  uses  BSMTH  to  calculate  each  of  the 
two  sub-splines.  Hence,  BSMCRV  calculates  a  smooth,  parametric  spline.  The  smoothing 
parameter  for  the  calls  to  BSMTH  is  determined  by  the  function  PRERR  and  the  stiffness 
weights  are  determined  by  the  subroutine  WTIBEG.  In  addition,  the  arc-length  is 
normalized  by  the  total  length  of  the  curve:  that  is,  the  parameter  used  is  not  the  arc- 
length  but  the  fractional  arc  length  along  the  curve.  Thus  the  parameter  s  varies 
between  0  and  1. 

An  example  of  a  spline  generated  by  BSMCRV  is  shown  in  Figure  1 1 .  Although  the 
data  points  show  a  large  amount  of  scatter,  an  excellent,  smooth  curve  has  been  found  to 
fit  the  data.  Notice  that  the  crossing  of  the  curve  over  itself  is  of  no  consequence  to 
BSMCRV. 


6  CONCLUDING  REMARKS 

The  computer  subroutines  presented  in  this  memorandum  extend  the  available 
libraries  of  spline  subroutines  at  DREA.  The  versatility  of  BSMTH  in  comparison  with  the 
BSPLIN  subroutine  SMOOTH  and  the  IMSL  subroutine  ICSSCU  ,  make  it  suitable  for  use  with 
a  far  greater  variety  of  data  sets.  In  particular,  the  ability  to  choose  the  splihe  order, 
the  ability  to  vary  the  spline  knots  independent  of  the  data  points,  and  the  ability  to 
change  the  'stiffness'  of  the  spline  at  specific  locations  via  the  stiffness  weights,  6n, 
allow  the  user  far  greater  control  over  the  spline  than  is  possible  with  SMOOTH  or 
ICSSCU.  Nor  need  the  choice  of  inputs  for  BSMTH  be  overly  difficult.  The  subroutines 
PRERR,  WTIBEG  and  NEWWTI  allow  the  user  to  generate  reasonable  sets  of  default  values 
for  the  smoothing  factor,  S,  and  the  stiffness  weights,  8n,  input  to  BSMTH.  Finally,  the 
restriction  that  the  data  points  be  splined  by  a  function  is  relaxed  if  one  chooses  to  use 
the  subroutine  BSMCRV.  Thus,  the  subroutine  library  provides  a  smoothing  spline  which 
provides,  at  once,  both  ease  of  use  and  great  freedom  and  flexibility. 


4  DEFAULT  VALUES  FOR  THE  STIFFNESS  WEIGHTS 


As  with  the  smoothing  parameter,  it  is  often  not  convenient  for  the  user  to  input  the 
values  for  the  stiffness  weights,  Sn,  n  =  1,N-k+1.  Two  subroutines  are  provided  which 
calculate  reasonable  values  for  the  parameters.  The  first  subroutine,  WTIBEG,  uses  the 
data  points  to  calculate  the  8n.  The  second,  NEWWTI,  uses  the  spline  coefficients  of  a 
previously  spline  approximation  of  the  data  to  calculate  new  values  for  5n. 

Both  subroutines  use  the  same  principle.  Default  values  for  the  5n  are  chosen  by 
setting  Sn  equal  to  a  predicted  value  for 


d^x)] 

l  dx*"  J 


2 

dx 


The  contributions  from  each  knot  interval  to  the  functional  F*  are  then  nearly  equal  and 
the  smoothing  will  not  be  dominated  by  one  short  segment  of  the  curve. 

In  WTIBEG,  it  is  assumed  that  m  =  2.  The  second  derivative  of  the  spline  in  any  knot 
interval  may  then  be  approximated  by  the  second  partial  difference  between  data  points 
near  that  knot  interval.  That  is,  if  xj.1  <  x  <  xj+1  and  t„  <  x  <  t^,  then 


f*'(x) 


Vi 


- 


j-i 


Vj-1  -  Vj 

Xj+1  *  Xj 


yj  -  yj-i 

Xj  -  Xj., 


(4.1) 


NEWWTI  uses  a  previous  spline  approximation  of  the  data  to  approximate  the 
integral  of  the  m**  spline  derivatives  In  any  knot  interval.  The  mtt  derivative  of  the  spline 
is  calculated  at  each  of  the  knots  and  the  integral  approximated  from  the  linear 
interpolation  of  these  values.  This  yields  the  formula 

f  [f(m)(x)]2dx  *  HtWUoXI^CWl2  ♦  f<m)(W)f<m)(Wi>  ♦ 

'Wi 

[f(m)(Wi>32>  (4-2) 

If  the  k  s  m+2,  this  method  is  exact  since  the  m,h  derivative  of  the  spline  is  then  linear 
between  the  knots. 

A  demonstration  of  the  ability  of  WTIBEG  to  choose  appropriate  choces  for  the 
stiffness  weights  is  shown  by  the  comparison  of  the  splines  in  Figures  1  and  2.  As 
explained  in  Section  2.6,  the  only  effective  difference  in  the  calculations  of  these  two 
splines  is  the  variation  in  the  stiffness  weights. 
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If  f  is  suitably  smooth,  then  the  first  term  on  the  right  side  of  equation  (3.8)  remains 
small  as  m  increases,  while  the  second  term  increases  rapidly.  Thus,  for  sufficiently  large 
m  and  N, 

N  N  N  N 

D&0  Vj  *  2  *  2  DS°  <£i>  =  2  D»?  eF2  s  9nm)  0-9) 

j=1  j=1  j=1  j=1 

so  that  an  estimate  for  O'2  is 


N-m  /  N  \2 

n=1  Vj=1  / 


(3.10) 


is  easily  calculated  from 


D<rvE  d^Vj 


(3.11) 


That  is,  D^ej  is  the  mth  divided  difference  of  the  data  set  {0,0 . ej,...0,0}. 


Thus,  in  order  to  estimate  a2,  and  hence  S,  It  Is  only  necessary  to  have  a  method 
for  determining  a  sufficiently  large  m.  The  domination  of  the  divided  differences  by  the 

N-m  N-m 

errors  is  characterized  by  a  large  number  of  changes  in  sign  between  yj  and 

j=1  j=1 

yj-  ,f  dominated  by  the  errors,  these  values  should  be  distribute  randomly  so  that, 
on  average,  one  expects  (N-m)/2  sign  changes.  Smooth  data  should  have  far  fewer.  The 
number  of  sign  changes  in  the  divided  differences  is  therefore  used  as  a  criterion  for 
determining  when  the  error  is  dominant. 

Figures  6,  8  and  10  demonstrate  the  ability  of  PRERR  to  calculate  appropriate 
smoothing  parameters.  Figure  8  shows  a  data  set  obtained  from  measurement?  of  the 
variation  of  sound  speed  with  depth  in  the  Atlantic  Ocean,  as  splined  by  an  ordinary  cubic 
spline  (the  subroutine  CUBSPL  from  the  BSPLIN  library  was  used).  Figure  9  shows  the 
same  spline  with  the  data  points  removed  so  that  the  curve  may  be  seen  more  easily.  It 
can  be  seen  that  the  curve  is  not  smooth,  especially  near  x»  15.  Figure  1 0  shows  the 
same  data  splined  using  BSMTH  with  the  smoothing  parameter  calculated  by  PRERR.  The 
fit  to  the  data  is  still  excellent  but  the  spline  is  now  smooth. 
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If  the  relative  magnitudes  of  the  en  accurately  reflect  the  errors  of  the  data  collection 
process,  then  averaged  over  a  large  number  of  data  sets  the  average  values  of  each  wn 
will  be  equal. 

<wn>  =  cr  for  all  n  (3.3) 


Here  angle  brackets  denote  averaging  over  an  ensemble  of  similar  data  sets. 


Since  f(x)  is  assumed  to  be  a  smooth,  well-behaved  curve,  it  should  be  possible  to 
fit  a  spline  curve  to  It  with  high  accuracy.  Hence,  the  X2  of  the  “best"  spline  is 


(3.4) 


if  It  may  be  assumed  that  the  errors  in  the  data  points  are  uncorrelated  and  that  Np  is 
sufficiently  large. 

Let  {gr  j=1,N}  be  any  set  of  numbers.  The  m1*  divided  difference  of  {gj}  is  a  linear 
transformation  of  the  g;  defined  iteratively  by 

■  8nj  (3.5) 


l(m). 


pf  Xn+m"*n 


n  =  1  ,N-m 


(3.6) 


where  $nj  Is  the  Kronecker  delta.  By  the  Mean  Value  Theorem,  if  f(x)  Is  a  Cm  function, 
then  for  any  {xj,  J*1,N}  there  is  a  $  in  (x^x,^)  such  that 


j*1 


D(m) 


f(Xj) 


m! 


(3.7) 


Thus,  from  equation  (3.1)  one  obtains 


f(m)«) 


*  2  “SS 

J«i 


m! 


(3.8) 
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order  of  the  spline  is  4  and  the  smoothing  exponent  m  is  2.  Hence,  the  only  difference 
between  this  spline  and  the  spline  calculated  by  SMOOTH  arises  from  the  effect  of  the 
stiffness  weights.  These  weights  have  been  decreased  near  x  =  8  and  x  =  15  to  allow 
the  spline  to  bend  rapidly  there.  The  stiffness  has  also  been  increased  between  x  =  10 
and  x  =  12  to  flatten  the  top  of  the  curve.  Notice  the  absence  of  wiggles.  These 
stiffness  weights  were  determined  by  the  subroutine  WTIBEG  (See  Section  4  and 
Appendix  A.5). 

For  the  spline  shown  in  Figure  3,  the  order  of  the  spline  was  increased  to  6.  In 
Figures  4  and  5,  the  second  derivative  of  the  SMOOTH  spline  and  the  sixth  order  BSMTH 
spline  are  shown,  respectively.  Since  the  SMOOTH  spline  is  necessarily  of  fourth  order, 
its  second  derivative  is  piecewise  linear. 

The  spline  shown  in  Figure  6  was  obtained  by  decreasing  the  spline  order  to  3,  and 
reducing  the  knots  as  shown.  At  the  positions  of  the  double  knots,  the  spline  need  no 
longer  have  continuous  derivative.  Splines  with  discontinuous  derivatives  cannot  be 
obtained  front  SMOOTH  or  ICSSCU.  For  certain  data  sets  they  are  necessary  to  obtain  an 
accurate  fit:  for  example,  when  splining  a  ship  hull  with  a  chine.  The  spline  shown  in 
Figure  7  carries  this  idea  one  step  further.  At  the  triple  knots,  the  spline  is  no  longer 
continuous  at  all.  While  a  use  for  a  completely  discontinuous  spline  may  not  be  evident, 
this  example  does  serve  to  illustrate  the  versatility  of  the  subroutine  BSMTH. 


3  CALCULATION  OF  INPUT  VALUES  FOR  THE  SPLINE  X2 

The  smoothness  of  the  splines  determined  by  BSMTH,  SMOOTH  and  ICSSCU  is 
regulated  by  the  input  parameter  S,  the  value  of  the  X2  of  the  resulting  spline.  It  is  often 
not  convenient  for  the  user  to  supply  this  input  parameter,  nor  is  an  appropriate  value 
likely  to  be  known.  In  this  section  an  algorithm  is  described  which  yields  an  appropriate 
value  for  the  parameter  S,  given  the  set  of  data  points  to  be  splined  and  their  associated 
errors  and  assuming  that  the  errors  are  uncorrelated.  Since  statistical  methods  are  used, 
the  algorithm  works  best  when  there  are  more  than  1 5  data  points.  The  algorithm  is 
implemented  in  the  function  subroutine  PRERR. 

Let  (xn,yn),  n  *  1  ,Np  be  the  data  points  and  e„  their  associated  errors.  It  is  assumed 
that  the  data  may  be  derived  from  some  unknown  "smooth"  curve,  f(x),  so  that 

yn*<W  +  «n  (3-D 

is  the  actual  error  of  the  0th  data  point.  This  must  not  be  confused  with  en,  which  is  the 
error  of  the  0th  data  point  estimated  by  the  collector  of  the  data.  The  en  are  known  -,  the 
£„  are  not. 

The  actual  errors  En  may  be  expressed 
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2)  S  <  X2,  :  p1  is  too  high.  Therefore,  set  phi  =  1 ,  X2W  =  X2,  and 
p2  =  pmln.  After  X22  is  determined  there  are,  again,  two  possibilities: 

i)  S  >  X22  :  p2  is  too  low.  Set  P|0  =  Pmin  and  X2|0  =  X22.  p3  is  now 
determined  such  that  (p3,S)  lies  on  the  straight  line 
interpolating  (pl0,  X2|0)  and  (phI>X2h,). 

ii)  S  <  X22  :  p2  is  too  high.  However,  p  cannot  be  decreased 
below  pm|n.  Therefore,  the  iteration  terminates. 

b)  Once  (P|O,X20),  and  (phl,X2h,)  have  been  determined  the  iteration  proceeds 
as  follows: 

1)  If  |S  -  X2n|  <  S/10,  the  Iteration  terminates. 

2)  If  x2n  -  x2lo  >  X2W  -  X2,,  then  p^,  is  determined  such  that  (p^.S) 
lies  on  the  straight  line  interpolating  (pn,X2n)  and  (ph|IX2h,). 

3)  If  X2n  -  X2to  <  X2hl  -  X2n,  then  p^,  is  determined  such  that  (Pp+i.S) 
lies  on  the  straight  line  interpolating  (pl0,X2|0)  and  (p„,X2n). 

This  procedure,  though  somewhat  more  complicated  than  the  simple  secant 
procedure  used,  for  example,  in  the  BSPLIN  subroutine  SMOOTH  (see  Reference  1, 
chapter  14),  converges  much  more  rapidly. 


2.6  Examples  of  splines  calculated  by  BSMTH 

As  examples  of  the  versatility  of  BSMTH  in  comparison  with  the  BSPLIN  subroutine 
SMOOTH  (the  IMSL  subroutine  ICSSCU  gives  splines  very  similar  to  SMOOTH),  a  simple  set 
of  data  points  has  been  splined  using  both  SMOOTH  and  BSMTH.  The  input  values  for  the 

data  point  errors,  e(  and  the  spline  X2  was  the  same  in  all  cases.  These  inputs  completely 
determine  the  spline  calculated  by  SMOOTH.  However,  the  versatility  of  BSMTH  becomes 
apparent  when  one  examines  the  many  qualitatively  different  curves  which  can  be  made 
to  fit  the  data  using  BSMTH.  These  curves  are  plotted  in  Figures  1  to  7. 

Figure  1  shows  the  spline  calculated  by  SMOOTH.  Notice  the  wiggles  caused  by  the 
inability  of  the  spline  to  bend  rapidly  near  the  points  x  •  8  and  x«  15.  The  small  crosses 
below  the  curve  indicte  the  positions  of  the  breakpoints  or  knots  of  the  spline.  For 
SMOOTH,  these  are  necessarily  at  the  data  point  abscissae,  with  the  exception  of  the 
second  and  next  to  last  data  point. 

Figure  2  demonstrates  the  effect  of  the  stiffness  weights  in  BSMTH.  The  knots  for 
this  spline  were  placed  at  the  data  points  (as  are  the  breakpoints  used  by  SMOOTH).  The 
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v*  and  Y*2  are  calculated  in  the  subroutine  SETUPR  during  the  calculation  of  R,j. 
During  the  iteration  for  p,  the  X2  is  evaluated  using  equation  (2.17)  in  the  subroutine 


XSQC. 


2.5  The  iteration  for  p 

The  spline  calculated  by  BSMTH  is  required  to  have  a  X2  equal  to  S,  a  value  input  by 
the  user.  This  is  implemented  by  iterating  over  the  value  of  a  in  equation  (2.7)  until 
|X2  -  S|  <  S/10.  In  practice,  BSMTH  iterates  over  p,  defined  in  equation  (2.13)  rather 
than  a. 

As  a  increases  from  0  to  1 ,  the  X2  of  the  spline  minimizing  G*  increases  from  some 
minimum  value  to  some  maximum  value.  However,  although  the  linear  system  of  equation 
(2.12)  is  theoretically  invertible  for  any  a  in  (0,1),  Rjj  is  not  invertible,  and  ,  depending  on 
the  positions  of  the  knots  with  respect  to  the  data  points  (see  Reference  1 ,  chapter  1 3), 
Pij  might  not  be  invertible  either.  Hence,  as  a  approaches  0  or  1 ,  there  will  be  numerical 
difficulties  in  the  inversion  of  equation  (2.12).  For  this  reason,  the  allowed  range  of  a, 
and  therefore  p  is  restricted.  The  upper  and  lower  limits  for  p  are  denoted  pm|n  and  p^, 
respectively.  pm)n  is  given  the  default  value  of  0.001  and  p^  the  default  value  of  1 000. 
These  values  have  been  found  adequate  to  circumvent  any  numerical  difficulties  when 
using  BSMTH,  though  they  may  be  changed  If  desired. 

The  iteration  for  p  is  divided  into  two  steps. 

a)  First,  values  of  p  and  their  corresponding  X2ls  are  determined.  These  are 
denoted  (P|0,X2l0),  and  (pw,X2hl).  They  are  determined  as  follows. 

Let  pn  denote  the  n*  value  of  p  determined  and  X2n  the  corresponding 

X2.  The  initial  guess  for  p  is  p1  *  1.  The  linear  system  of  equation  (2.12) 

is  inverted,  and  the  X2  of  the  spline  Is  evaluated.  There  are  two 
possibilities: 

1)  S  >  X2,  :  In  this  case,  p,  is  too  low.  Set  plo  =  1  and  X2le  »  X2V  p2  is 
set  to  A  gain  there  are  two  cases: 

I)  S  >  X22  :  p2  is  still  too  low.  However,  p  cannot  be  increased 
above  pmax.  Therefore,  the  iteration  terminates. 

II)  S  <  X22  :  p2  is  too  high.  Set  pw  »  p2  and  X2hl  «  X22.  p3  is  now 
determined  such  that  (p3,S)  lies  on  the  straight  line 
Interpolating  (pte,  X2le)  and  (pw,X2h,). 
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2.3  Evaluation  of  P„j,  Rnj,  and  vn 

The  matrix  Pnj  is  evaluated  in  the  subroutine  SETUPP.  Since  B^fx)  is  a  piecewise 
polynomial  of  order  k,  the  integrals  in  the  definlton  of  Pnj  can  be  evaluated  by  a  series  of 
integrations  by  parts. 

np  k-m 

pnj  =  2  M  (2.1  4) 

p=1  q=1 


If  k  >  2m,  then  m-q  will  become  negative.  By  convention  B$P  (x),  for  q  <  0,  is  defined  to 
be  the  qth  integral  of  Bjk(x).  The  subroutine  BSPLVD,  from  the  BSPLIN  library,  is  used  to 
evaluate  the  derivatives  of  the  B-splines.  If  k  >  2m,  integrals  of  the  B-spiines  must  also 
be  calculated.  This  is  most  easily  accomplished  by  calculating  the  coefficients  of  the 
knot  sequence  corresponding  to  the  integral  of  each  B-spline  (see  Reference  1,  page 
1 50)  and  then  using  the  BSPLIN  subroutine  BVALUE  to  evaluate  It.  However,  to  calculate 
the  spline  coefficients,  k-2m  knots  must  be  appended  to  the  knot  sequence.  Thus  the 
dimension  of  the  array  containing  the  knots  is  required  to  be  Nk  +  max(0,k-2m). 

Owing  to  the  left  continuity  of  the  B-splines  as  implemented  in  the  subroutines 
BSPLVD  and  BVALUE,  and  the  discontinuity  of  the  higher  derivatives  of  the  B-splines,  they 
cannot  be  evaluated  right  at  the  knots.  Instead,  they  are  evaluated  at 
(0.9999t|  +  .0001 1)+1)  and  (0.000 1  tj  +  .9999t,+1)  for  each  knot  interval. 

In  practice,  P,j  in  equation  (2.12)  is  replaced  by  Py/A,  where  A  is  a  normalizing 
factor  used  to  ensure  that  the  elements  of  P,j  are  of  order  1.  This  averts  unwanted 

overflows  and  underflows.  It  has  no  effect  on  the  minimization  of  6*  as  the  factor  A  can 
be  absorbed  into  a  redefinition  of  a.  A  is  defined  by 


A  « 


(N-k+D2”1*1 


(2.15) 


The  matrix  R,j  is  evaluated  in  the  subroutine  SETUPR  making  use  of  the  BSPLIN 
library  subroutine  BSPLVB  to  evaluate  Bltk(xn).  This  subroutine  is  a  modification  of  the 
subroutine  L2APPR  In  the  BSPLIN  library. 
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for  given  a,  and  iterating  over  a  until  the  spline  has  the  required  X2. 


2.2  Minimization  of  G*  for  given  a 

Using  equations  (2.1),  (2.2),  (2.5),  and  (2.6),  the  functional  G*  may  be  written  in  the 
following  form: 

N  N  N 

((1  -a)Rnj  ♦  aPnj)^j  -  2(1-p)£vA  +  (1-p)Y2  (2.7) 

n=1 j=1  n=1 


where 


R 


pj 


"p 

£ 


Bp.k(xn)Bj,k(xn) 


(2.8) 


,  .  B  V  6p  f  ^  BP.k(xn)  dm 

PJ  &  dx" 


.Bj.kCXn) 


dx 


ynBj.k(Xn) 


Y2 


(2.9) 


(2.10) 


(2.11) 


G*  Is  minimized  with  respect  to  the  spline  coefficients,  0n,  when 
N 

£  (Rni  *  Ppni#i  8  Vn  (2.12) 

J-1 

where 


P  ■ 


a 

1  -a 


(2.13) 


Since  both  Py  and  Ry  are  symmetric,  banded,  positive  definite  matrices,  the  subroutines 
BCHFAC  and  BCHSLV  In  the  BSPLIN  library  are  appropriate  for  the  solution  of  the  linear 
system  in  equation  (2.12). 
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X2  measures  the  degree  to  which  the  curve  approximates  the  data  and  is  minimized  when 
the  curve  interpolates  the  data. 

The  second  functional  is  a  measure  of  the  smoothness  of  the  spline  function.  Its 
definition  relies  on  the  observation  that,  for  a  smooth  function,  the  average  values  of  its 
high  order  derivatives  will  be  considerably  lower  than  those  of  a  'wiggly'  function.  Hence, 
one  uses  the  functional 


(2.3) 


as  a  measure  of  the  smoothness  of  the  spline. 


The  spline  desired  is  that  which  has  a  given  X2  while  minimizing  F.  In  practice,  this 
is  found  by  finding  the  spline  which  minimizes  the  functional 

G  =  aX2  +  (1  -a)F  (2.4) 

for  given  a.  Since  G  is  quadratic  in  the  spline  coefficients  |9n,  this  amounts  to  the  solution 

of  a  linear  system.  An  iteration  is  then  done  to  find  the  value  of  a  for  which  X2  has  the 
required  value.  As  implemented  by  Reinsch  and  de  Boor,  the  splines  are  necessarily 
cubic,  and  the  knots  are  constrained  to  be  the  data  point  abscissae,  x„,  n=1,N. 

A  shortcoming  of  the  above  algorithm  Is  that  the  second  derivative  of  the  spline  is 
minimized  even  in  places  where  one  might  expect  It  to  be  high:  that  is,  where  the  data 
shows  a  pronounced  bend.  This  problem  has  been  avoided  in  BSMTH  by  generalizing  the 
functional  F  to 


(2.5) 


The  basis  for  the  linear  space  of  spline  functions  has  been  chosen  to  be  the  B-spline 
basis  (see  Reference  1 ,  chapter  8).  The  n^  B-spline  of  order  k  is  denoted  Bn  k(x)  and  the 
knots  of  the  spline  are  denoted  t„,  n  *  1  ,Nk. 

The  coefficients  £„  can  be  used  to  alter  the  'stiffness'  of  the  spline  between  the 
pair  of  knots,(tB_w.1,tn_h+2).  In  regions  where  the  spline  curve  is  required  to  be  very  flat, 
the  £n  will  be  large.  In  regions  where  the  spline  is  expected  to  have  high  curvature,  the 
£„  will  be  small. 


The  required  spline  is  found  by  minimizing  the  functional 
G*  ■  aX2  ♦  (1-a)F* 


(2.6) 


2 


Section  1 


decreased.  Two  subroutines,  WTIBEG  and  WTINEW  (see  Section  4)  are 
provided  which  calculate  appropriate  default  values  for  the  stiffness 
weights  from  the  data  points  or  from  previous  spline  fits  to  the  data, 
respectively. 

5)  SMOOTH  and  ICSSCU  implement  smoothing  by  minimizing  the  second 
derivative  of  the  the  spline.  BSMTH  allows  one  to  choose  the  derivative 
which  is  to  be  minimized,  again  allowing  more  control  over  the  character  of 
the  spline. 

BSMTH  does  have  the  drawback  that  it  is  somewhat  slower  than  the  other 
subroutines,  though  usually  at  most  by  a  factor  of  two.  However,  much  of  the  extra  time 
can  often  be  made  up  by  reducing  the  number  of  knots  of  the  spline  with  no  deterioration 
in  the  quality  of  the  fit  (the  execution  time  is  roughly  proportional  to  the  number  of 
knots).  Moreover,  when  splining  in  two  dimensions,  the  time  savings  involved  in  having 
the  data  points  independent  of  the  knots  far  outweigh  the  slight  inefficiency  of  BSMTH. 

In  the  following  sections  the  algorithms  for  each  of  the  subroutines  is  discussed  in 
detail.  User's  guides  including  sample  runs  of  the  subroutines  are  given  in  Appendix  A. 
The  computer  code  for  each  subroutine  is  given  in  Appendix  B. 


2  THE  BSMTH  ALGORITHM 


2.1  Implementation  of  Smoothing 


The  technique  used  in  ICSSCU  and  SMOOTH,  for  constructing  a  smooth  spline  curve 
through  a  given  set  of  data  is  an  extension  of  an  algorithm  first  proposed  by  Whittaker4 
and  later  considered  by  Schoenberg5,  Reinsch3  and  de  Boor1.  The  idea  is  to  define  two 
functionals  dependent  quadratically  on  the  spline  coefficients  for  which  one  is  solving. 
One  functional  is  the  X2  of  the  spline  curve, 


(2.1) 


where 

(Xn.yn),  n*1  ,Np;  are  the  data  points  to  be  interpolated, 
en  is  the  error  associated  with  the  n-th  data  point,  and 
N 

f(x)  «  0nf„(x)  (2.2) 

n*1 


The  functions  fn(x)  are  the  basis  functions  for  the  linear  space  of  spline  functions.  The 
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=  5,  If  0.9*X2  of  the  spline  of  the  data  point  abscissae  is  greater  than 
the  value  predicted  by  PRERR. 

=  6,  If  1.1  *X2  of  the  spline  of  the  data  point  abscissae  is  less  than 
the  value  predicted  by  PRERR. 

=  7,  8,  9  ,  As  for  IER  =  4,  5,  and  6  respectively,  but  for  the  spline  of 
the  data  point  ordinates. 

BCOEFX :  An  array  of  length  N  containing  the  B-spline  coefficients  of  the 
spline  of  the  abscissae. 

BCOEFY :  An  array  of  length  N  containing  the  B-spline  coefficients  of  the 
spline  of  the  ordinates. 

ARCL  :  An  array  of  length  N  containing  the  estimated  arc-length  at  each 
data  point. 

Via  COMMON  /  CHISQ  / 

XSQ  =  X2  of  the  spline 
WORK  SPACE 

WTI  :  An  array  of  length  N  which  is  used  to  contain  the  stiffness  weights 
as  calculated  by  WTI8EG. 

G  :  An  array  of  length  NPT*IMAX  used  as  work  space  in  the  function 
PRERR. 

R  :  An  array  of  length  3*K*N 

IWK  =  max(NKT1  ,K**2) 

WK  :  An  array  of  length  4*IWK 
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The  following  data  has  been  splined  using  BSMCRV.  The  resulting  spline  has  been 
plotted  in  Figure  1 1 . 


NPT  - 

22,  N  - 

20,  K  -  4, 

NKT1  - 

J 

X  ( J) 

Y  ( J) 

E  (J) 

1 

-0.24 

-0.09 

0.3 

2 

0.06 

0.03 

0.3 

3 

0.20 

0.18 

0.3 

4 

0.31 

0.39 

0.3 

5 

0.43 

0.47 

0.3 

6 

0.45 

0.63 

0.3 

7 

0.39 

0.77 

0.3 

8 

0.38 

0.86 

0.3 

9 

0.29 

0.94 

0.3 

10 

0.11 

0.97 

0.3 

11 

0.06 

1.03 

0.3 

12 

-0.08 

1.02 

0.3 

13 

-0.17 

1.00 

0.3 

14 

-0.23 

0.97 

0.3 

15 

-0.31 

0.91 

0.3 

16 

-0.29 

0.80 

0.3 

17 

-0.33 

0.70 

0.3 

18 

-0.30 

0.59 

0.3 

19 

-0.15 

0.45 

0.3 

20 

0.06 

0.33 

0.3 

21 

0.26 

0.11 

0.3 

22 

0.66 

0.03 

0.3 

The  spline  coefficients  and  the  fractional  arc  length  values  returned  by 
BSMCRV  are 
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J 

BCOEFX(J) 

BCOEFY(J) 

ARCL(J) 

1 

-0.35219 

-0.20888 

0 . 00000E+00 

2 

-0.19B60 

-0.10975 

0. 90250E-01 

3 

-0. 39770E-01 

-0. 1964GE-01 

0.14756 

4 

0.10954 

0.11428 

0.21378 

5 

0.24B99 

0.27045 

0.25406 

B 

0.37S81 

0.43389 

0.29910 

7 

0.4B3B3 

0.61444 

0.34165 

8 

0.40S11 

0.79B60 

0.36694 

9 

0.24539 

0.97128 

0.40057 

10 

0. 5B733E-01 

1.0224 

0.45154 

11 

-0. 13699 

1 . 0425 

0.47336 

12 

-0.32258 

0.88117 

0.51256 

13 

-0. 34446 

0.B986B 

0.53832 

14 

-0.21200 

0.532S8 

0.55705 

15 

-0.53392E-01 

0.39341 

0.58498 

IB 

0.11140 

0.26571 

0.61621 

17 

0.27779 

0.14252 

0.64630 

18 

0.44BB9 

0. 43413E-01 

0.67814 

19 

0.62444 

0. 14056E-01 

0.73545 

20 

21 

22 

0.8000B 

-0.3181SE-01 

0.80301 

0.88606 

1 . 0008 

The  values  returned  via  COMMON  /  CHISQ  /  are 

XSQX  «  0. 7G300E-01 ,  XSQY  -  0.11847,  SX  -  0.71920E-01,  SY  . 


0.117BB 
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A.2  BSMTH  :  User's  Guide 

SUBROUTINE  BSMTH(S,JDER,NPT,X,Y,E,N,K,NKT1  J.WTI.BCOEF.R.IWK.WK.IER) 

PURPOSE:  BSMTH  calculates  the  spline  of  order  K,  with  knots  T(I),I=1,NKT  which  has  chi- 

square  of  S  with  respect  to  the  data  points  X(I),Y(I),I=1,NPT,  and  which  has  as  small 

a  JDER,h  derivative  as  possible. 

LANGUAGE:  FORTRAN 

USAGE:  EXECUTE  mainpgm,BSPLIN:HLLYSP/LIB,BSPLIN:BSPLIN/LIB 

CALLS  subroutines  SETUPQ,  SETUPR,  XSQC,  SMODAV  and  INTERV,  BCHFAC,  BCHSLV  from 

the  BSPLIN  library 

INPUT 

S  :  The  chi-square  of  the  spline  with  respect  to  the  data  will  be  within 
10%  of  S,  if  possible.  As  S  is  increased  the  spline  becomes 
smoother  but  farther  from  the  data  points.  Function  PRERR  can  be 
used  to  give  a  value  for  S  if  a  reasonable  value  is  not  known. 

JDER  :  The  integral  of  the  square  of  the  JDER*1  derivative  of  the  spline  is 
minimized  (subject  to  the  constraint  that  XSQ  =  S).  If  smooth  curves 
are  desired  a  value  of  JOER  *  2  is  appropriate.  JDER  should  be  non¬ 
negative  and  less  than  K. 

NPT  :  The  number  of  data  points. 

X  :  An  array  of  length  NPT  containing  the  data  point  abscissae  in 
ascending  order. 

Y  :  An  array  of  length  NPT  containing  the  data  point  ordinates. 

E  :  The  errors  of  the  data  points.  The  smaller  the  error  the  closer  the 

spline  will  come  to  that  point. 

N  :  The  number  of  B-spllnes  used  to  represent  the  spline. 

K  :  The  order  of  the  spline. 

NKT1  *  N  +  K  +  max(0,K-2*JDER) 

T  :  An  array  of  length  NKT1  the  first  N+K  elements  of  which  contain  the 
knot  sequence  (in  ascending  order).  The  remaining  array  elements 
are  used  In  subroutine  SETUPP. 

WTI  :  An  array  of  length  N  of  which  only  the  first  N-K+1  elements  are  used 
(rather  than  passing  In  an  otherwise  superfluous  argument).  WTI(I) 
is  a  weight  for  the  integral  of  the  square  of  the  JOER*1  derivative  of 
the  spline  between  T(l+K-1)  and  T(l+K).  The  larger  WTI(I)  is  the 
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smoother  the  integral  will  be  over  this  region.  These  weights  are 
relative:  i.e.  changing  all  the  WTI  by  a  constant  factor  will  not 
affect  the  resulting  spline. 

IER  =  0,  If  JDER,T,WTI  and  the  first  N#K  elements  of  R  are  as  on  the 
previous  call  to  BSMTH  (  this  means  that  the  matrix  P  need  not  be 
recalculated  ) 

=  1 ,  if  P  is  to  be  recalculated 
Via  COMMON  /  PLIMS  / 

PMIN  =  Minimum  allowed  value  of  p  (See  (Section  2.5)).  Default  is  1.0E-03 
PMAX  =  Maximum  allowed  value  of  p.  Default  is  1.0E+03. 

OUTPUT 

IER  *  0,  Calculation  has  been  successful 

*  1 ,  If  JDER  >  K  -  1 

*  2,  If  NKT1  <  N  +  K  +  max(0,K-2*JDER) 

=  3,  If  IWK  <  max(NKT1  ,K**2) 

*  4,  If  more  than  30  Iterations  are  required  to  find  the  correct  value 

for  P.  Indicates  numerical  difficulties  in  the  solution  of  the  linear 
system 

=  5,  If  the  X2  of  the  spline  >  1 .1*S 

*  6,  If  the  X2  of  the  spline  <  .9*S 

BCOEF  :  An  array  of  length  N  containing  the  B-spline  coefficients  of  the 
spline. 

Via  COMMON  /  CHISQ  / 

XSQ  *  X2  of  the  spline 
WORK  SPACE 

R  :  An  array  of  length  3»K*N 

IWK  «  max(NKT1  ,K**2) 

WK  ;  An  array  of  length  4*IWK 
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The  following  input  data  has  been  splined  using  BSMTH.  A  plot  of  the  spline  is 
shown  in  Figure  2. 


-  10.0 

.  JDER 

-  2  ,  NPT 

-  1 0.X 

-  4  ,  N 

-  10  ,  NKT  -  . 

J 

X(J) 

Y(J) 

E(J) 

T(J) 

UTI (J) 

1 

7.0 

0.0 

0.005 

7.0 

0. 51183E-02 

2 

8.0 

13.5 

0.5 

7.0 

0.55789 

3 

9.0 

15.5 

0.5 

7.0 

1.1778 

4 

10.0 

14.5 

0.5 

7.0 

1.1778 

5 

11.0 

15.5 

0.5 

9.0 

2.6500 

6 

12.0 

15.0 

0.5 

10.0 

0.88333 

7 

13.0 

14.5 

0.5 

11.0 

0. 15407E-01 

8 

14.0 

15.0 

0.5 

12.0 

9 

15.0 

13.5 

0.5 

13.0 

10 

16.0 

0.0 

0.005 

14.0 

11 

16.1 

12 

16.1 

13 

16.1 

14 

16.1 

The  spline  coefficients  obtained  were 


J  BCOEF  (J) 

1  0. 3195009E-04 

2  1 A. 03531 

3  14.93398 

4  15. 00996 

5  15.11922 

6  15.09243 

7  15.00040 

8  14.90065 

9  13.10809 

10  -2.075629 

and  the  X2  of  the  spline  was 


XSQ  -  9.9912 


~Tr 
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A.3  NEWWTI  :  User's  Guide 

SUBROUTINE  NEWWTI(NOLD,BCOEF,NKTOLD, TOLD, NKTNEW.TNEW, NWTI, WTI.JDER) 

PURPOSE:  NEWWTI  uses  a  previously  calculated  spline  fit  to  predict  values  for  the 
stiffness  weights  5,  for  use  in  BSMTH. 

LANGUAGE:  FORTRAN 

USAGE:  EXECUTE  mainpgm,BSPLIN:HLLYSP/UB,  BSPLIN:BSPLIN/LIB 
CALLS  subroutines  SMODAV  and  BVALUE  from  the  BSPLIN  library. 

INPUT 

NOLD  :  Number  of  B-splines  for  old  spline  fit. 

BCOEF  :  Array  of  length  N  containing  the  B-spline  coefficients  for  the  old 
spline  fit. 

NKTOLD  :  Number  of  knots  for  the  old  spline  fit. 

TOLD  :  Array  of  length  NKT  containing  the  knots  for  the  Old  spline  fit. 

NKTNEW:  Number  of  knots  for  the  new  spline  fit. 

TNEW  :  Array  of  length  NKT  containing  the  knots  for  the  new  spline  fit. 

NWTI  :  Number  of  stiffness  weights  for  the  new  spline  fit. 

JDER  :  The  order  of  derivative  minimized  by  BSMTH. 

OUTPUT 

WTI  :  Array  of  length  NWTI  containing  the  stiffness  weights,  5,. 
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The  following  input  data  has  been  used  to  generate  stiffness  weights  by  NEWWT. 
This  data  is  the  output  data  from  the  example  in  Section  A.2. 

NOLD  -  4  ,  NKTOLD  -  14  .  NKTNEU  -  14  ,  NUTI  -  7  ,  JDER  -  2 

J  TOLD  (J)  TNEU(J) 


1 

7.0 

7.0 

2 

7.0 

7.0 

3 

7.0 

7.0 

4 

7.0 

7.0 

5 

9.0 

9.0 

G 

10.0 

10.0 

7 

11.0 

11.0 

8 

12.0 

12.0 

3 

13.0 

13.0 

10 

14.0 

14.0 

11 

1G.1 

IE.  1 

12 

IB.  1 

IB.  1 

13 

1G.1 

IB.  1 

14 

IB.  1 

16. 1 

The  stiffness  weights  obtained  were 

J  UTI (J) 

1  0. 10000E-02 

2  0.5G836E-01 

3  1.159G 

4  0.51772 

5  4.8525 

E  0. 14507E-01 

7  0.10000E-02 
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A.4  PRERR  :  User's  Guide 

FUNCTION  PRERR(NPT,X,Y,E,IMAX,G,WK,IFLAG) 

PURPOSE:  This  function  calculates  the  mean  error  in  the  data  points.  The  smoothing 
parameter  used  by  BSMTH  may  then  be  determined  by:  S  =  NPT#PRERR**2. 

LANGUAGE: FORTRAN 

USAGE:  EXECUTE  mainpgm,BSPLIN:HLLYSP/LIB 

CALLS  subroutines  PARDIF 

INPUT 

NPT  :  The  number  of  data  points. 

X  .  An  array  of  length  NPT  containing  the  data  point  abscissae  in 
ascending  order. 

Y  :  An  array  of  length  NPT  containing  the  data  point  ordinates. 

E  :  The  errors  of  the  data  points. 

IMAX  :  The  maximum  number  of  partial  differences  taken  is  2*!MAX.  The 

suggested  value  for  IMAX  is  5. 

IFLAG  =  0,  If  the  calculation  is  to  be  done  from  scratch. 

-  1 ,  If  X,  E  and  G  have  not  been  changed  since  the  previous  call. 

OUTPUT 

PRERR  returns  the  mean  error  in  the  data  points. 

WORK  SPACE 

WK  :  An  array  of  length  NPT 

G  :  An  array  of  length  N*IMAX 
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The  mean  error  in  the  following  data  has  been  predicted  by  PRERR.  Splines  of  this 
data  are  shown  in  Figure  8,  9,  and  10  and  are  discussed  in  Section  3. 

N  -  40,  1I1AX  -  5 


J 

X(J) 

V(J) 

J 

X(J) 

YU) 

1 

0.00 

1507.89 

21 

22.70 

1482.57 

2 

3.G3 

1507.85 

22 

23.00 

1481.83 

3 

7.2G 

1507.81 

23 

23.50 

1480.34 

4 

10.90 

1507.77 

24 

24.20 

1478.83 

5 

12.20 

1507.17 

25 

2G.10 

1476.17 

6 

13.70 

1505.02 

2G 

27.00 

1475.02 

7 

14.10 

1502.49 

27 

28.10 

1472.68 

8 

14.50 

1501.53 

28 

29.00 

1470.30 

9 

14.80 

1499.50 

29 

29.90 

1468.70 

10 

15.20 

1498.27 

30 

30.  G0 

1467.92 

11 

15.30 

149G.94 

31 

44.00 

1463.25 

12 

15.40 

1495.93 

32 

52.70 

1464.54 

13 

15.70 

1494.92 

33 

58.40 

1466.09 

14 

1G.G0 

1492.87 

34 

65.20 

1469.74 

15 

1G.80 

1491.84 

35 

74.50 

1475.43 

1G 

17.30 

1490.09 

36 

80.30 

1478.08 

17 

18.40 

1488.33 

37 

94.50 

1482.65 

18 

20.90 

148G.22 

38 

110.00 

1487.18 

19 

21.80 

1484.77 

39 

119.10 

1489.40 

20 

22.50 

1483.31 

40 

158.10 

1491.40 

PRERR  returned  the  value  0.20361 
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A.S  WTIBEG  :  User's  Guide 

SUBROUTINE  WTIBEG(NPT,X,Y,NKT  J.NWTI.WTI) 

PURPOSE:  WTIBEG  uses  the  data  points  to  calculate  values  for  the  stiffness  weights  6, 
for  use  in  BSMTH. 

LANGUAGE:  FORTRAN 

USAGE:  EXECUTE  mainpgm,BSPLIN:HLLYSP/LlB,  BSPLIN:BSPLIN/LIB 
CALLS  subroutines  SMODAV  and  BVALUE  from  the  BSPLIN  library. 

INPUT 

NPT  :  The  number  of  data  points. 

X  :  An  array  of  length  NPT  containing  the  data  point  abscissae  in 
ascending  order. 

Y  :  An  array  of  length  NPT  containing  the  data  point  ordinates. 

NKT  :  Number  of  knots  for  the  spline. 

T  :  Array  of  length  NKT  containing  the  knots  for  the  spline. 

NWTI  :  Number  of  stiffness  weights  for  the  spline  fit.  NWTI  =  NKT-2*K+1 

where  K  is  the  order  of  the  spline. 

OUTPUT 

WTI  :  Array  of  length  NWTI  containing  the  stiffness  weights,  5,. 
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The  following  input  data  has  been  used  to  generate  stiffness  weights  in  WTIBE6. 
The  spline  obtained  from  this  data  is  discussed  in  Section  2,6  and  is  plotted  in  Figure  2. 


NPT  =10,  NKT  =14,  NWTI  =  7 


J 

X(J) 

V(J) 

T  (J) 

1 

7.0 

0.0 

7.0 

2 

8.0 

13.5 

7.0 

3 

9.0 

15.5 

7.0 

4 

10.0 

14.5 

7.0 

5 

11.0 

15.5 

9.0 

6 

12.0 

15.0 

10.0 

7 

13.0 

14.5 

11.0 

8 

14.0 

15.0 

12.0 

9 

15.0 

13.5 

13.0 

10 

16.0 

0.0 

14.0 

11 

16.1 

12 

16.1 

13 

16.1 

14 

16.1 

The  stiffness  weights  obtained  were 


a 

1  0. 51183E-02 

2  0.55789 

3  1.1778 

4  1.1778 

5  2.6500 

6  0.88333 

7  0. 15407E-01 


nnoooooonooorjoooooonoooooorjoc) 
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Appendix  B 
SUBROUTINE  LISTINGS 


C  )|c 

C  *  THESE  COnPUTER  SUBROUTINES  ARE  THE  PROPERTY  OF  THE  * 

C  *  CANADIAN  DEPARTMENT  OF  NATIONAL  OEFENCE  _  * 

C  9|c  $ 

C  *  THEY  SHALL  BE  USED  ONLY  FOR  PURPOSES  AUTHORISED  * 

C  *  BY  THE  DEPARTMENT  .  * 

C  £ 

C  *  THEY  SHALL  NOT  BE  DISCLOSED  TO  A  THIRD  PARTY  * 

C  *  UITHOUT  THE  WRITTEN  PERMISSION  OF  THE  * 

C  *  DEPARTMENT  .  * 

C  *  # 


B.1  BSMCRV 


SUBROUT I NE  BSMCRV (NPT , X , Y , E , N , K . NKT1 , T , UT I , BCOEFX , BCOEFY , 

*  R ,  I  UK ,  UK ,  ARCL ,  G ,  I ER ) 

- C 

c 

Given  data  points  (X(I), Y(I)),  I«1,NPT  BSMCRV  finds  a  smooth  C 

curve  approximating  them  by  sp I i n i ng  the  abscissae  and  ordinates  C 
separately  with  respect  to  the  fractional  arc- length  along  the  C 
spline.  An  approximation  for  the  arc  length  at  each  point  is  C 

obtained  from  the  distances  between  the  points.  BSMTH  is  used  to  C 
spline  the  absissae  and  the  ordinates.  PRERR  is  used  to  C 

determine  a  smoothing  factor  for  the  splines  and  UTIBEG  is  used  C 
to  determine  stiffness  weights.  C 

C 

AUTHOR:  David  Hally  ,  May  1381  C 

C 

USAGE:  C 

EXECUTE  ma i npgm , BSPL I N : HLLYSP/L I B , BSPL I N : BSPL I N/L I B  C 

C 

CALLS  PRERR, BSMTH, UTIBEG  C 

C 

INPUT:  C 

C 

VIA  SUBROUTINE  ARGUMENTS:  C 

C 

NPT  :  The  no.  of  data  points  C 

X  :  An  array  of  length  NPT  containing  the  data  point  C 

Y  :  An  array  of  length  NPT  containing  the  data  point  C 

ordinates.  C 

E  :  The  errors  of  the  data  points.  The  smaller  C 

the  error  the  closer  the  spline  will  come  C 

to  that  point.  C 
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N  :  The  no.  of  B-splines.  C 

K  :  Tne  order  of  the  spline.  C 

NKTl  «  N-t-K-t-max (0, K-2*JDER)  C 

T  :  An  array  of  length  NKTl  the  first  N+K  elements  of  C 

which  contain  the  knot  sequence.  The  variable  used  C 
to  parametrize  the  curve  is  the  arc  length  divided  by  C 
the  total  length  of  the  curve.  Thus  the  knots  must  C 
span  the  interval  [0,11.  Default  gives  a  uniform  C 

distribution  of  knots  over  this  interval.  C 

C 

IER  -  0  If  defaults  are  desired  C 

»  1  If  defaults  are  not  desired  C 

C 

COMMON  /NODFLT/  C 

C 

IMAX  :  2*IMAX  is  the  max.  no.  of  divided  differences  allowed  C 

to  find  the  error  (used  in  function  PRERR) .  Default  C 

value  is  5  C 

SMFACT  :  See  comments  belou  C 

C 

COMMON  /INTEXP/  :  C 

C 

JDER  :  The  integral  of  the  square  of  the  JDER-th  derivative  C 

of  the  spline  is  minimized  (  subject  to  the  CON-  C 

straint  that  XSQ-S) .  Default  value  is  2  .  C 

JOER  must  not  exceed  K-l  C 

DEFAULTS:  C 

C 

If  IER  «  0  on  input  then:  C 

JOER  -  2  C 

SMFACT  -  1.0  C 

IMAX  -  5  C 

T ( I )  «  (I -K) / (N-K+l) ,  I-l.NKT  i.e.  knots  are  uniformly  C 

distributed  in  (0,1)  C 

C 

OUTPUT:  C 

c 

IER  -  0  .  Iteration  converged  C 

-  1  ,  If  JDER  >  K-l  C 

-  2  .  If  NKTl  <  N+K  +MAX ( 0 , K -2*JDER )  C 

-  3  ,  If  I UK  <  max (NKTl, K**2)  C 

■  4  ,  If  iteration  for  PI  in  BSMTH  did  not  converge  C 

during  the  spline  of  the  X-values  C 

-  5  ,  If  the  chi-square  of  the  spline  of  the  X-values  C 

returned  by  BSMTH  >  1.1*S  (i.e.  PMAX  in  BSMTH  C 

i s  too  sma II)  C 

■  S  ,  If  the  chi-square  of  the  spline  of  the  X-values  C 

returned  by  BSMTH  <  .9*5  (i.e.  PMIN  in  BSMTH  C 

is  too  large)  C 

«  7,8,3  As  for  IER-4,5, S, respect ive I y,  but  for  the  spline  C 
of  the  Y-values  C 

BCOEFX  »  Array  of  length  N  containing  the  B-spline  coefs.  for  C 
the  X-values  of  the  curve  C 
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B.6  SETUPP 

SUBROUT  I NE  SETUPP (NPT , E , JDER , T, NKT , N, K , UT I , P , A , DB , DB1 , UK , A1 ) 


SETUPP  calculates  the  matrix  P  (  see  Ref.  Manual  ) 
AUTHOR:  David  Hally  ,  May.  1S81 
CALLED  by  BSMTH 
CALLS  SMOOAV 

from  BSPLIN  library  BVALUE, BSPLVD 


REAL  T  (NKT) ,  P  (K,  N) ,  A  (K,  K)  ,DB  (K,  K) ,  DB1  (K,K)  ,UT1  (N),  UK  (NKT) , 

*  E(NPT) ,A1 (K,K) ,T1,T2,H,HI 

INTEGER  NKT, JDER, N.K.MMAX, I , J.L.LJl, IKJ.M, II. KM 

C  A  normalizing  factor  H  is  calculated.  Normalization  by  H  ensures  that 
C  most  of  the  elements  of  P  are  of  order  1. 

H-  ( (T  (N+K)  -T  (in  /FLOAT  (N-K+l) )  ** (2*JDER-1 ) 

C  P  is  initialized  and  extra  points  are  added  to  the  knot  sequence 
C  to  allow  the  calculation  of  higher  order  B-sp lines  if  necessary 
C  in  the  integration  by  parts. 

H-H/ (SMOOAV (NPT, E) **2*SM0DAV (N-K+l , UTI ) ) 

00  10  I-i.K 
00  10  J-l.N 
10  P ( I , J) “0. 

IF (N+K. EQ. NKT) GO  TO  30 
DO  20  J-NKT, N+K+l , -1 
20  T(J)-T (N+K) *1.0001 

30  MMAX-M I N0 ( JDER , K -JDER ) 

C  An  iteration  over  the  intervals  between  knots  is  begun. 

00  140  I-K.N 

IF(T(I+1) .EQ. T ( I ) ) GO  TO  140 
Il-I-K+1 

C  The  derivatives  of  the  B-splines  needed  in  the  integration  by  parts 
C  are  calculated  using  BSPLVD.  Due  to  the  left  continuity  of  BSPLVD 
C  the  derivatives  are  evaluated  close  to  but  not  right  at  the  knots. 

T1«.9999*T (I )+. 0001 *T (1+1) 

T2- . 0001*T ( 1 ) +. 9999*T U  +1 ) 

CALL  BSPLVD (T.K, T1 , I , A.DB.K) 

CALL  BSPLVD (T,K, T2, I , A.0B1 , K) 


C  The  integrals  of  the  B-splines  needed  in  the  integration  by  parts 


UO(JCJL)L)LHJDL)L)L) 
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CALL  PARDIF(N,X.UK,  J.  J+l.N-J) 

J«  J+l 
J2-J2+2 

IF(NSGNCH.GT.D/2. ) GO  TO  130 
IF (J.GE. 1  MAX— 2) GO  TO  120 
GO  TO  70 

120  SDEV. (D/2. -NSGNCH) /SORT (D) 

C  The  error  is  determined  from  the  divided  difference  by  taking  the 
C  root  mean  square  of  the  divided  difference  values  weighted  by 
C  the  expected  value  for  a  unit  error  (given  by  GU.J+l)). 

C  anomalously  high  values  are  discarded  and  the  resulting  error 
C  is  corrected  by  multiplying  prerr  by  1.14 

130  DEV-0.0 

DO  140  I -J+l.N-J 

140  DEV- (UK (I ) /G ( 1 , J+l ) ) **2+DEV 

I F  (DEV.  EQ. 0.0) RETURN 
NM2J-N-J2 

PRERR-SQRT  (DEV/FLOAT  (NH2J)  )*2.0 
DO  150  I -J+l.N-J 

IF  (ABS  (UK  (I )  /GU.J+l) )  .LT.  PRERR)  GO  TO  150 
DEV-DEV- (UK ( I ) /G ( I , J+l ) ) **2 
NI12J-Nn2J-l 
150  CONTINUE 

PRERR-SQRT  (DEV/FLOAT  (NFI2J) )  *1 . 14 

RETURN 

END 


Appendix  6 


41 


REAL  X  (N) ,  Y  (N)  ,E  (N)  ,UK  (N)  ,G  (N, I  MAX) ,DEV,D,PRERR 
INTEGER  NSGNCH , NN2  J ,  N ,  K , I HAX ,  J ,  J2 ,  KI1 1 N , KflAX , I  FLAG, I 

COfinON  /CERR/  SDEV 

INAX-niN0(N/2,inAX) 

SDEV.0.0 

PRERR-0.0 

IF ( I  FLAG. EQ. 1 ) GO  TO  50 


C  G  is  calculated. 

DO  10  I-l.N 
G(I.l)-EU) 

DO  10  J-2,  II1AX 
10  G  (I ,  J) -0. 0 

DO  30  I-l.N 
DO  20  J-l.N 

20  UK (J) -0.0 

UK(l)-Ed) 

DO  30  J-l.IMAX-l 

KniN-MAX0(I-J-l,J) 

KnAX-niN0(I+J+l,N-J+l) 

CALL  PARD I F  (N ,  X ,  UK ,  J-l ,  KP1 1 N ,  KtlAX) 

IF  (I  .GT.  J)Kf1IN-KniN+l 
IF  d+J.LT.N+1) KMAX-KMAX-l 
DO  30  K-KniN,KnAX 

30  G  (K,  J+l)  -UK  (K)  **2+G  (K,  J+l) 

DO  40  J-l .  iriAX-1 
DO  40  I-J+l.N-J 

40  G ( I , J+l ) -SQRT (G ( I , J+l ) ) 

C  Divided  differences  are  taken  until  the  no.  of  sign  changes  is 
C  greater  than  that  expected  for  random  data.  If  II1AX-2  iterations 
C  occur  first  SDEV  is  set  to  the  number  of  standard  deviations 
C  that  NSGNCH  is  belou  its  expected  value. 

50  DO  G0  I-l.N 
60  UK  (I) -Yd) 

J-0 

J2-0 


C  The  no.  of  sign  changes  in  the  divided  differences  is  determined. 

70  NSGNCH-0 

I -J+l 

80  DO  90  K-I+l.N-J 

IF  (UK  (I  )*UK  (K) )  100,90, 110 
90  CONTINUE 

100  NSGNCH-NSGNCH+1 

110  I-K 

IF  (I ,LT.N-J)G0  TO  80 

D-N-J2-1 
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B.S  PRERR 

FUNCTION  PRERR (N.X.Y.E. I  MAX, G, UK, I  FLAG) 


This  subroutine  calculates  the  mean  error  in  the  data  points  (X,Y) 
by  taking  divided  differences  until  the  no.  of  sign  changes  in 
the  I-th  divided  difference  is  that  expected  from  random  data. 

The  error  is  then  determined  by  assuming  that  the  contribution 
from  the  smooth  curve  underlying  the  data  is  negligible. 

AUTHOR:  David  Hally  .  Jan.  1981 

USAGE: 

EXECUTE  ma i npgm . BSPL I N : HLL  YSP/L I B 
CALLS  PARDIF 
INPUT  : 

N  »  No.  of  data  points 

X  :  An  array  of  length  N  containing  the  data  point 

abscissae  in  ascending  order. 

Y  :  An  array  of  length  N  containing  the  data  point 

ordinates. 

E  :  An  array  of  length  N  containing  the  relative  errors  of 

the  data  points.  The  absolute  errors  are  obtained  by 
multiplying  the  returned  value  of  PRERR  by  the 
relative  errors. 

I  FLAG  «  0  ,  If  calculation  is  to  be  done  from  scratch 

-  1  ,  If  If1AX,X,E,  and  G  have  the  same  value  as  in  the 

previous  cal  I 

G  -  Array  of  dimensions  N.IMAX.  G(I,J)  is  the  expectation 

value  of  the  J-th  divided  difference  given  an  error 
of  E(I)  in  the  I-th  data' point.  If  IER-0  G  Is 
calculated  ;  otherwise  it  is  assumed  known. 

IMAX  :  2*IMAX  is  the  max.  no.  of  divided  differences  alloued 
IMAX  -  5  is  suggested 

OUTPUT  : 

PRERR  -  The  calculated  mean  error  in  the  data 
VIA  COTTON  /  CERR  / 

SDEV  :  The  no.  of  sign  changes  in  the  divided  difference  used 
to  calculate  PRERR  is  greater  than  that  expected  for 
random  data  less  SDEV  standard  deviations 

WORK  SPACE  : 

UK ( I )  OF  DIMENSION  N 


U(JU(JUL)UU(JUL)U(JUUUUL)(JU(J(JUUL)L)UCJL)L)UU(JL)UUL)UUUUUL)(J(JUU(JL)UU 
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B.4  PARDIF 

SUBROUTINE  PARDIF (N, X, F, J, IrtlN, I  MAX) 


PARDIF  calculates  the  divided  difference  cf  the  data  pcints 
(X  (I ) ,  F  (I ) ) ,  I  —  IMIN,  IflAX.  Tc  avcid  ever-  or  underflows  the  X 
intervals  are  normalized  by  the  factor  H- (X (N) -X  (1) ) /N  .  This  is 
of  no  consequence  in  PRERR  since  only  ratios  of  partial  diff¬ 
erences  are  of  significance. 

AUTHOR:  David  Hally  ,  Hay.  1981 

CALLED  by  PRERR 


REAL  X (N) ,F (N) , H 
INTEGER  J.iniN.iriAX.I.N.IT 

H-  (X(N)-X(l))  /FLOAT  (N) 

DO  10  IT-1,2 

DO  10  1-IflIN, INAX-IT 

10  F  (I )  -H*(F  ( 1+1 )  -F  (I ) )  /  (X  (I+I T)  -X(I-J) ) 

DO  20  I -IMAX-2, IMIN, -1 
20  F ( I +1 ) -F ( I ) 

f  tiniN)  -0.0 
F  UMAX) -0.0 
RETURN 
END 
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B2-BV ALUE (TOLD , BCOEF , NOLD , KOLD , T2 , JDER) 

UTI (IU) =  (Bl#  (B1+B2)  +B2**2)  *  (TNEU  ( I +KNEU)  -TNEU  (I+KNEU-1) )  /3. 
10  CONTINUE 

C  The  modal  average  of  LIT  I  is  determined  and  LIT  1(1)  is  set  to 
C  UTIAV/UTI (I ) 

UTIAV-SMODAV(NUTl.UTl) 

IF  (UTIAV.EQ.0.0)GO  TO  30 
UMIN-UTI AV#1 . 0E-03 
UMAX«UT  I  AV#1 . 0E+03 
DO  20  IU-l.NUTI 
dummy-uti ( i w) 

IF  ( (LIT  I  ( I  LI) .  GT.  Util  N)  .AND.  (UTI  (ID)  . LT. UMAX) ) UTI  (IU). 

*  UTIAV/UTI (IU) 

IF (DUMMY. GT. UMAX) UTI (IU) -1 . 0E-03 
20  IF (DUMMY. LE.UMINIUTI (IU)-1.0E+03 

RETURN 

30  DO  40  IU-l.NUTI 
40  UTI <IU)-1.0 

RETURN 
END 
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B.3  NEWWTI 

SUBROUTINE  NEUUTI (NOLD. BCOEF, NKTOLD. TOLD. NKTNEU, TNEU, NUTI . UTI , 
*  JOER) 


NEUUTI  uses  the  previous  spline  fit  to  predict  values  for  the 
integral  weights  UTI  for  use  in  BST1TH. 

AUTHOR:  David  Hally  ,  Aug.  1381 

USAGE: 

EXECUTE  mai npgm, BSPL I N : HLL  YSP/L I B 


CALLS  SnODAV 

from  BSPLIN  library  :  BVALUE 

INPUT: 


NOLD 

BCOEF 

NKTOLD 

TOLD 

NKTNEU 

TNEU 

NUTI 

JDER 

OUTPUT: 


»  No.  of  B-sp lines  for  old  spline  fit 
«  Array  of  length  N  containing  the  B-spline  coefficients 
for  the  old  spline  fit 
■  No.  of  knots  for  the  old  spline  fit 

-  Array  of  length  NKT  containing  the  knots  for  the  old 

sp line  fit 

-  No.  of  knots  for  the  new  spline  fit 

-  Array  of  length  NKT  containing  the  knots  for  the  new 

spline  fit 

-  No.  of  integral  weights  for  the  new  spline  fit 
•  The  order  of  derivative  minimized  by  BSHTH 


UTI  -  Array  of  length  NUTI  containing  the  integral  weights 


REAL  BCOEF (NOLD) .TOLD (NKTOLD) , TNEU (NKTNEU) , UTI (NUTI) , 
*  Bl,B2,Tl,T2.UniN,UnAX,UTIAV,DUnMY 

I NTEGER  NOLD . NKTOLD , NKTNEU , KOLO , KNEU .NUTI, JDER , I , I U 

K0LD-NKT0L0-N0LD 
KNEU- (NKTNEU-NUT I +1 ) /2 
IU-0 


C  On  each  knot  interval  the  integral  of  the  square  of  the  JDER-th 
C  derivative  of  the  given  spline  is  approximated 

DO  10  I -1. NUTI 

IF  (TNEU  (I+KNEU-1)  .EQ.  TNEUd+KNEU)  )G0  TO  10 
IU-IU+1 

T1 - . 9999*TNEU ( I +KNEU-1 ) +. 0001 *TNEU ( I +KNEU) 
T2».0001*TNEU ( I +KNEU-1 ) +. 9999*TNEU (I+KNEU) 

B1 -BVALUE ( TOLD . BCOEF , NOLD . KDLD , T 1 , JDER ) 
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PHI -PI 

P1-.1*PL0+.9*PHI 
GO  TO  220 

C  Similarly,  i f  PI  is  very  close  to  PLO,  it  is  possible  that  XSCkXSQLO 
C  In  this  case  PLO  is  set  to  PI,  XSQLO  to  XSQ  and  PI  to  . 9*P1+,1#PH1 

170  IF (XSQ. GT. XSQLO) GO  TO  180 

XSQLO-XSQ 
PL0-P1 

P1-.9*PL0+.1#PHI 
GO  TO  220 

180  IF ( (S-XSQLO) ,LT, (XSQH1-S) )G0  TO  190 

P2- (PI -PHI) * (S-XSQHI ) / (XSQ-XSQHI ) +PHI 
GO  TO  200 

190  P2«  (PI -PLO) # (S-XSQLO ) / (XSQ-XSQLO) +PLD 

200  IF  (P2.LT.PDG0  TO  210 

PLO-P1 
XSQLO-XSQ 
P1-P2 

IF (Pl.GT.PHI )P1» (PLO+PH I) 11 . 

GO  TO  220 

210  PHI -PI 

XSQHI-XSQ 

P1-P2 

IF (PI .LT. PLO) PI- (PLO+PH I ) /2. 

220  CONTINUE 

IER-4 


C  BCOEF  is  returned  to  its  correct  value  (see  comment  before  call  to 
C  XSQC) . 

230  DO  240  I-l.N 

240  BCOEF ( I ) -BCOEF ( I ) +UK (1,1) 

RETURN 

END 
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90  BCOEF ( I ) -BCOEF ( I ) -UK (1,1) 

XSQ-XSQC (N,  K, BCOEF ,R  (1, 1 ,2)  ,UK  (1,3)  ,UK  (1,4) , YSQ) 

C  If  XSQ  is  with  in  .  1*S  of  S  the  iteration  terminates. 

C  The  first  value  of  XSQ  calculated  is  for  s  Pl-1.  ,  then  Pl-PMIN  or 
C  Pl-PHAX  depending  on  whether  S  is  less  or  greater  than  XSQ.  The  third 
C  value  of  PI  is  predicted  by  linear  interpolation  of  the  tuo  known 
C  points.  The  known  P’s  and  their  correspond! ng  XSQ’s  are  then: 

C  (PLO, XSQLO) , (PHI , XSQHI ) ,  and  (PI, XSQ)  respectively.  Subsequently 
C  improved  values  of  PI  are  predicted  by  a  linear  interpolation 
C  of  (PI, XSQ)  and  either  (PLO, XSQLO)  or  (PHI, XSQHI)  depending  on 
C  whether  S  is  closer  to  XSQLO  or  XSQHI. 

C  If  XSQHI  <  S  or  XSQLO  >  S  initially  the  iteration  terminates. 

100  '  IF(ABS(S-XSQ).LT.S*.1)G0  TO  230 

GO  TO (110, 130), IT 
GO  TO  1S0 

110  IF (S.LT.XSQ)GO  TO  120 

XSQLO-XSQ 
PL0-P1 
Pl-PMAX 
GO  TO  220 

120  XSQHI -XSQ 

PHI -PI 
PI -PI1  IN 
GO  TO  220 

130  IF (Pl.EQ.PniN)GO  TO  140 

IF (S.LE.XSQ)GO  TO  135 
IER-5 
GO  TO  230 

135  XSQHI -XSQ 

PHI -PI 
GO  TO  150 

140  IF (S.GE.XSQ)GO  TO  145 

IER-S 
GO  TO  230 

145  XSQLO-XSQ 

PL0-P1 

150  PI - (PHI -PLO) * (S-XSQLO) / (XSQHI -XSQLO) +PL0 

GO  TO  220 

C  It  is  possible  that  due  to  numerical  inaccuracy  in  the  evaluation 
C  of  XSQ,  that  XSQ>XSQHI.  This  would  normally  only  occur  if  PI  iB 

C  very  close  to  PHI.  Hence  PHI  is  set  to  PI,  XSQHI  to  XSQ  and 

C  PI  to  . 1*PL0+. 9*PHI 


160 


IF (XSQ. LT. XSQHI )G0  TO  170 
XSQHI -XSQ 
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IER-1 

RETURN 

10  IF  (NKTl.GE.NKT+f1AX0(0,K-2*JOER)  )G0  TO  20 
IER-2 
RETURN 

20  IF ( (IUK.GE.NKT1) .AND. (IUK.GE. K#*2) ) GO  TO  30 
IER-3 
RETURN 


C  The  matrices  P  and  R  are  calculated  in 
C  SETUPP  and  SETUPR  respectively. 

30  IF (1ER.EQ.0)GO  TO  40 

CALL  SETUPP (NPT, E. JOER, T, NKT1, N.K.UTI.RQ, 1,3), UK. UK (1.2). 
*  UK  (1,3)  ,UK  (1,4)  ,R) 


C  The  array  UK(.,1)  is  determined  so  that  UK(.,1)  approximates  Y. 

C  This  is  necessary  for  accurate  calculation  of  XSQ. 

40  IER-0 

DO  B0  I-l.N-1 
DYSQ-0. 

DO  50  J-l.K-1 

50  DYSQ-DYSQ+T (I+J) 

DYSQ-DYSQ/FLOAT (K-l) 

CALL  INTERV (X, NPT, DYSQ, LEFT, I1FLAG) 

I F  (flFLAG.  EQ.  1 )  LEFT  -NPT -1 

DYSQ- (DYSQ-X (LEFT) ) / (X (lEFT+1 ) -X (LEFT) ) 

80  UK ( I . 1 ) -Y (LEFT) * (1 . -DYSQ) +Y (LEFT+1 ) *DYSQ 

UK (N,l)-Y (NPT) 

CALL  SETUPR (NKT, T,N,K,NPT, X, Y,E, Y5Q,R (1, 1 ,2) ,UK,UK (1,2)  ,UK (1,3) ) 


C  An  iteration  is  begun  which  changes  PI  until  XSQ  is  within 
C  .1*S  of  S 

Pl-1. 

XSQLO-0.0 
XSQHI-8.0 
DO  220  IT-1,30 
DO  70  I-1,K 
DO  70  J-l , N-I+l 

70  R (I , J, 1) »P1*R(I , J,3)+R(I , J,2) 

C  The  equation  R»BCOEF-VCT  is  solved  by  first  finding  the 
C  Choleeky  factorization  of  R,  then  by  solving  for  BCOEF. 

CALL  BCHFAC(R,K,N,UK (1,4) ) 

DO  80  I-1,N 

80  BCOEF ( I )-UK (I ,2) 

CALL  BCHSLV(R,K,N, BCOEF) 


C  The  chi-square  of  the  solution  is  determined 
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over  this  region.  These  ueights  are  relative:  i.e. 
changing  all  the  UTI  by  a  constant  factor  ui I  I  not 
affect  the  resulting  spline. 

IER  »  0  ,  If  JDER, T, UTI  and  the  first  N*K  elements  of  R  are 
as  on  the  previous  call  to  BSMTH  (  this  means 
that  the  matrix  P  need  not  be  recalculated  ) 

«  1  ,  if  P  is  to  be  recalculated 

VIA  COnnON  /  PLIMS  /  : 

PAIN  «  Min.  allowed  value  of  PI  (  See  comment  describing 

iteration  for  correct  chi-square  ). Default  is  1.0E-03 
PMAX  «  Max.  allowed  value  of  PI.  Default  is  1.0E+03. 

OUTPUT: 


IER  «  0  ,  Calculation  has  been  successful 

-  1  ,  If  JDER  >  K-l 

-  2  .  If  NKT1  <  N+K+max (0, K-2*JDER) 

-  3  ,  If  I  UK  <  max(NKTl,K**2) 

-  4  ,  If  more  than  30  iterations  are  required  to  find  the 

correct  value  for  P.  Indicates  numerical  difficul¬ 
ties  in  the  solution  of  the  linear  system 

-  5  ,  If  the  chi-square  of  the  spline  >  1.1*S 

«  S  ,  If  the  chi-square  of  the  spline  <  .9*S 

BCOEF  :  An  array  of  length  N  containing  the  B-spline  coeffi¬ 
cients  of  the  spline. 

VIA  COMMON  /  CHISQ  /  : 


XSQ  -  the  chi-square  of  the  spline 
UORK  SPACE: 

R  :  An  array  of  length  3*K*N 

I UK  -  max (NKT1,K*#2) 

UK  :  An  array  of  length  4*IUK 


REAL  BCOEF (N).T(NKTl), UTI (N),R(K,N,3) ,UK(IUK,4), 

*  X  (NPT) ,  Y  (NPT) ,  E  (NPT) , 

*  PI , P2 , PH I . PLO, XSQ , XSQHI . XSQLO . YSO , ALF , DYSQ, S 

I NTEGER  N , K , NKT , NKT1 , JDER , NPT , MFLAG , LEFT , I  UK , I ER , I T , I . J 

COMMON  /  CHISQ  /  XSQ.DUM(3) 

COMMON  /  PLIMS  /  PM IN, PMAX 

DATA  PMIN  /  1.0E-03  /.PMAX  /  1.0E+03  / 

C  The  input  data  is  checked  for  simple  errors 


NKT-N+K 

IF (JDER.LT.K)GO  TO  10 
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B.2  BSMTH 

SUBROUT  I NE  BSMTH  (S ,  JDER , NPT , X , Y , E , N , K , NKT 1 , T , UT I ,  BCOEF , 
*  R, IUK.UK, IER) 


BSMTH  calculates  the  spline  of  order  K,  with  knots  T(I),I-1,NKT 
which  has  chi-square  of  S  with  respect  to  the  data  points 
X  (I ) ,  Y  (I ) ,  1  «1,NPT,  and  which  has  as  small  a  JDER-th  derivative 
as  possible. 

AUTHOR:  David  Hally  ,  May  1981 
USAGE: 

EXECUTE  ma i npgm , BSPL I N : HLLYSP/LIB . BSPL I N : BSPL IN/LIB 

CALLS  SETUPQ . SETUPR , XSQC 

from  BSPL I N  library:  I NTERV , BCHFAC , BCHSLV 


INPUT: 

S 


JDER 


NPT 

X 

Y 

E 


N 

< 

NKT1 

T 


UTI 


The  chi-square  of  the  spline  with  respect  to  the  data 
will  be  within  10%  of  S,  if  possible.  As  S  is 
increased  the  spline  becomes  smoother  but  farther 
from  the  data  points.  Function  PRERR  can  be  used 
to  give  a  value  for  S  if  a  reasonable  value  is  not 
known. 

The  integral  of  the  square  of  the  JDER-th  derivative 
of  the  spline  is  minimized  (  subject  to  the  con¬ 
straint  that  XSQ-S).  If  smooth  curves  are  desired 
a  value  of  JDER-2  is  appropriate.  JDER  should  be 
non-negative  and  less  than  K. 

The  no.  of  data  points 

An  array  of  length  NPT  containing  the  data  point 
abscissae  in  ascending  order. 

An  array  of  length  NPT  containing  the  data  point 
ordinates. 

The  errors  of  the  data  points.  The  smaller 
the  error  the  closer  the  spline  will  corns 
to  that  point. 

The  no.  of  B-sp lines. 

The  order  of  the  spline. 

N+K+max (0, K-2*JDER) 

An  array  of  length  NKT1  the  first  N+K  elements  of 

which  contain  the  knot  sequence  (in  ascending  order), 
the  remaining  array  elements  are  used  in  subroutine 
SETUPP. 

An  array  of  length  N  of  which  only  the  first  N-K+i 
elements  are  used  (rather  than  passing  in  an  other¬ 
wise  superfluous  argument).  UTI  (I)  is  a  weight 
for  the  integral  of  the  square  of  the  JDER-th  deriv¬ 
ative  of  the  spline  between  TU+K-l)  and  T(l+K).  The 
larger  WT 1(1)  is  the  smaller  the  integral  will  be 
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1  MAX-5 
JDER-2 

DO  40  1-1, N+K 

40  T  (I) -FLOAT  (I-K) /FLOAT  (NUTI) 

C  The  error  in  the  X-values  are  found  by  calling  the  function 
C  PRERR  and  the  integral  weights,  UT I,  by  calling  UTIBEO. 

C  They  are  splined  using  BSMTH  using  fractional  arc  length 
C  to  parametrize  the  data  points.  Similarly  for  the  Y-values. 

C  NOTE:  The  parameter  SM  to  be  used  in  BSMTH  should  be 
C  expected  to  be  NPT#PRERR#*2 . 

C  However,  due  to  the  sensi  tivi  tyof  parametric  splines  to 
C  data  error,  it  has  been  found  that  slightly  higher  values  of 
C  SM  sometimes  give  better  results.  SMFACT  has  been  included  as 
C  a  knob  to  increase  (or  decrease)  SM  :  SM«SMFACT*NPT*PRERR#*2  . 
C  default  value  for  SMFACT  is  1.0. 

150  IER-1 

CALL  UT I BEG (NPT , ARCL , X , N+K . T , NUT I , UT I ) 

SX«SMFACT*PRERR (NPT, ARCL, X.E, I  MAX, G, UK, 0>**2*FLDAT (NPT) 

CALL  BSMTH (SX , JDER , NPT , ARCL . X , E . N , K . NKT1 , T , UT 1 , BCDEFX , 

*  R, 1UK.UK, 1ER) 

XSQX-XSGY 

IF ( (IER.NE.0) .AND. (1ER.LT.4) ) RETURN 
IER-1 

S Y-SMF ACT*PRERR (NP  T , ARCL , Y , E , I MAX , G , UK , 1 ) **2*FL0AT (NPT) 

CALL  UT I  BEG (NPT , ARCL , Y , N+K , T , NUT I , UTI ) 

CALL  BSMTH (S Y , JDER , NPT , ARCL , Y , E , N , K , NKT1 , T , UT I , BCOEFY , 

*  R, IUK.UK, IER) 

RETURN 

END 
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BCOEFY  «  Array  of  length  N  containing  the  B-spline  coefs.  for 
the  Y-values  of  the  curve 

ARCL(I)  «  Arc  length  at  the  I-th  data  point/total  length  of  curve 

via  connoN  /  chisq  / 

XSQX  ■  Chi-square  of  the  spline  of  the  abscissae 
XSQY  -  Chi-square  of  the  spline  of  the  ordinates 
SX  -  Required  Chi-square  of  the  abscissae  (  as  determined  by 
PRERR  ) 

SY  -  Required  Chi-square  of  the  ordinates  (  as  determined  by 
PRERR  ) 

via  common  /crvlth/  : 

SNEUL  :  The  total  arc  length  of  the  curve 
UORK  SPACE  : 

R  j  An  array  of  length  3*K#N 

I UK  -  max (NKT1,K**2) 

UK  :  An  array  of  length  4*IUK 

G  :  An  array  of  dimensions  NPT, IMAX  (used  by  PRERR) 

UTI  i  Array  of  length  N  used  for  the  integral  weights  for 
BSHTH 


REAL  X  (NPT) ,  Y  (NPT) .  E  (NPT) ,  ARCL  (NIPT) ,  BCOEFY  (N) , BCOEFX  (N) , 

*  T (NKT1) , UTI (N) ,G(NPT, IMAX) ,UK(IUK,4) ,R(K,N,3) , 

*  SHF ACT 

INTEGER  NPT.N.K.NKTl ,NUTI , IUK, IER, IHAX, I ,K1 , IU, ID 

COMMON  /NODFLT/  SMFACT, I  MAX 
COMMON  /INTEXP/  JOER 
COMMON  /CHISQ/  XSQY, XSQX, SY.SX 

NUTI-N-K+1 

C  ARCL(I) , I-l.N  is  initialized  by  connecting  the  data  points  with 
C  straight  lines. 

ARCL  (1) -0.0 
DO  10  I -2, NPT 

10  ARCL  (I )  -ARCL  (I  -D+SQRT  ( (X  (I )  -X  ( I  -1) )  **2+  (Y  (I )  -Y  (I  -1)  )**2) 

OLDL-ARCL (NPT) 

DO  20  I -2, NPT 

20  ARCL ( I ) -ARCL ( I ) /OLDL 

C  If  IER.NE.0  non-default  values  of  SMFACT, IMAX  and  IMAX  are 
C  taken  from  the  COMMON  block  /NODFLT/  and  JDER  from  COMMON  /INTEXP/ 

IF (1ER.NE.0)GO  TO  150 
SMFACT -1.0 


uuuuuuuuuuuuuuuuuuuuuuuuuuu 


44 


Appendix  B 


C  are  calculated  by  calculating  the  coefs.  of  the  knot  sequence 
C  correspond i ng  to  the  integral  of  each  B-spline  and  then  calling 
C  BVALUE  to  evaluate  these  at  the  appropriate  points. 

1 F  12# JOER .  GE . K ) GO  TO  90 
DO  80  J-l.K 
IKJ-I-K+J 
UK(IKJ)-1. 

IF (J.EQ.K)GO  TO  50 
DO  40  L-IKJ+1, 1+1 
40  UK  (L)  -0.0 

50  DO  70  M«1,K-2*JDER 

KM-K+M-1 
00  50  L-IKJ, 1+1 

UK  (L)  -UK  (L)  *  (T  (L+KM)  -T  (L) )  /FLOAT  (KM) 

I F  (L.  NE .  1 )  UK  (L)  -UK  (L)  +UK  (L-l ) 

50  CONTINUE 

A1 (J,M)-BVALUE(T,UK.N.K+M,T2,0) 

70  A(J,M)-BVALUE(T,UK,N,K+M,T1,0) 

80  UK (IKJ) -0. 0 

C  The  elements  of  P  are  determined  by  integration  by  parts. 

90  DO  130  L-l.K 

DO  130  J-l.L 
LJ1-L-J+1 
IKJ-I-K+J 
HI  -H 

IF (MMAX.LT.l)GO  TO  110 
DO  100  M-l.MMAX 

P  (LJ1 . 1 K J)  -P (I _ II ,  IKJ)+HI#UTI  (ID* (DB1  (L .  JDER-M+1 ) * 

*  DB1 (J, JOER+M) -OB (L, JDER-M+1 )*DB(J, JDER+M) ) 

100  HI— HI 

110  IF (K.LE.2*JDER)G0  TO  130 

DO  120  M-1,K-2*JDER 

P(LJ1, IKJ)-P(LJ1, IKJ)+HI*UTI (I1)*(A1  (L,M)*0B1 (J. JDER 

*  +nnAX+M) -A (L, M) *DB ( J, JDER+MMAX+M) ) 

120  HI— HI 

130  CONTINUE 

140  CONTINUE 
RETURN 
END 


onoooonooooonr) 


Appendix  B 


45 


B.7  SETUPR 

SUBROUT I NE  SETUPR (NKT,T,N,K,NPT,X, Y,E, YSQ,q, Yl, VCT, VCT1) 


The  matrix  R, the  arrays  Yl.VCT  and  VCT1,  and  the  number  YSQ  are 
ca I cu I ated 

(  SETUPR  is  based  closely  on  L2APPR  by  Carl  de  Boor, in 
A  Practical  Guide  to  Splines,  p.  255) 

AUTHOR:  David  Hally  ,  May.  1381 

CALLED  BY  BSMTH 

CALLS  from  BSPLIN  library  BSPLVB 


REAL  T(NKT) ,R(K,N) ,VCT(N) ,BIATX(20) ,X(NPT) ,Y(NPT) ,E(NPT) ,DU, 
#  Yl (N) ,  VCT1 (N) , YSQ, DYSQ 

I NTEGER  N ,  K ,  NXT ,  NPT ,  LEFT ,  LEFTfIX ,  I ,  J ,  1111 ,  J  J ,  LL 
YSQ-0. 

DO  20  J-l.N 
VCT1  (J)-0. 

VCT(J)-0. 

DO  10  I-l.K 
R(I,J)-0. 

10  CONTINUE 

20  CONTINUE 

C  The  LL-th  data  point  is  positioned  within  the  knot  sequence. 

LEFT-K 
LEFTNX-0 
DO  80  LL-l.NPT 

30  IF(LEFT.EQ.N)GO  TO  40 

IF(X(LL) .LT.T (LEFT+1) ) GO  TO  40 

LEFT-LEFT+1 

LEFTHK-LEFTflK+l 

GO  TO  30 

C  R  is  calculated  by  calling  BSPLVB  to  evaluate  the  B-sp lines  at 
C  the  data  points. 

40  CALL  BSPLVB (T,K,1,X(LL) ,LEFT,BIATX) 

DYSQ-Y(LL) 

00  50  MII-l.K 

DYSQ-DYSQ-BIATX(nn)*Yl(LEFT-K+nn) 

50  CONTINUE 

do  70  nn-i,x 

DU-BI  ATX  (fin)  /E  (LL)  **2 
J-LEFTflK+nn 

VCT1  ( J) -VCT1 ( J) +0YSQ#0U 
VCT  (J)  »DU#Y  (LLJ+VCT  (J) 


>  T 
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1-1 

00  S3  JJ-MM,K 

R  ( I .  J)  -B I  ATX  ( JJ)  *0UI+R  ( I ,  J) 
I-I+l 


63 

CONTINUE 

70 

CONTINUE 

YSQ-  (DYSQ/E  (LL ) )  **2+YSQ 

80 

CONTINUE 

RETURN 

END 
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B.8  SMODAV 

FUNCTION  SNODAV(NPT.X) 


SNOOAV  returns  a  modal  average  of  the  numbers  in  X 
AUTHOR  :  David  Hally  ,  Aug.  1981 
USAGE  : 

EXECUTE  ma i n-pgm , BSPLIN: HLLYSP/L I B 
INPUT  : 

NPT  -  No.  of  values  to  be  averaged 

X  -  Array  of  length  NPT  containing  values  to  be  averaged 
RETURNS: 

SNODAV  -  Nodal  average  of  the  values  in  X 


REAL  X (NPT) , XBOX (11 ) . SUNBOX (10) , XNI N. XNAX. XRATI 0. SCALE. SNODAV 
INTEGER  I BOX (10) ,NPT, ISUH.NBOX, I ,  J 

C  The  range  of  the  values  is  found  and  broken  into  NBOX  logarithmic 
C  intervals,  such  that  the  ratio  of  the  smallest  to  the  largest 
C  possible  no.  in  each  interval  does  not  exceed  NPT,  but  also 
C  subject  to  the  constraint  2  <  NBOX  <  11. 

XNIN-1 . 0E+30 
XNAX-1.0E-30 
00  10  I-l.NPT 

IF (X (I ) .LE. 0.0)GO  TO  10 
XNAX-ANAX1  (X(I).XNAX) 

XNIN-ANIN1 (X(I).XNIN) 

10  CONTINUE 

I F (XN I N . EQ. 1 . 0E+30 ) GO  TO  90 
XRATI 0-XNAX/XN I N 

NBOX-ALOG10 (XRATI 0) /ALOG10 (FLOAT (NPT) ) 

NBOX-N 1 N0  (NBOX ,  1 8 ,  NP  T  /5 ) 

NBOX-NAX0 ( 3 , NBOX ) 

SC ALE-XR A  T 1 0** ( 1 . /NBOX ) 

XBOX (1) -XNI N 
XBOX (NBOX+1) -XNAX 


C  The  no.  of  X-values  uithin  each  interval  is  calculated 

00  28  I -1. NBOX 
SUNBOX (I ) -0.0 
20  I  BOX ( I ) — 0 

00  30  I -2, NBOX 

XBOX ( I ) -XBOX ( I -1 ) *SCALE 


30 
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DO  S3  Ul.NPT 

DO  40  J-2.NB0X+1 

40  IP (X ( I ) . LE. XBOX (J) ) GO  TO  50 

50  SUMBOX ( J-l ) -SUnBOX  (J-l 

50  I  BOX (J-l) -I  BOX (J-l )+l 

C  Denote  by  Xmid  the  X-va I ue  such  that  there  are  an  equal  no.  of  X-values 
C  Doth  smaller  and  greater  than  Xmid.  SHODAV  is  the  average  value  of  all 
C  the  X’s  in  the  interval  containing  Xmid. 

I sun-0 

DO  70  Ul.NBOX 

isun-isun+iBOXd) 

I F  ( I  Sun . GE . NP  T /2 ) GO  TO  80 
70  CONTINUE 

80  snooAv-sunBOxin/iBoxd) 

RETURN 


90 


SnODAV-0.0 

RETURN 

END 
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B.9  WTIBEG 

SUBROUT  I NE  UTIBEG (NPT, X, Y, NKT, T, NUTI , UTI ) 


UTIBEG  uses  the  data  points  to  predict  values  for  the  integral 
ueights  UITI  for  use  in  BSI1TH. 

AUTHOR:  David  Hally  ,  Aug.  1981 

USAGE: 

EXECUTE  mai npgm, BSPL I N : HLL YSP/L I B 


INPUT: 

NPT 

X 

Y 

NKT 

T 

NUTI 


OUTPUT: 

UITI 


-  No.  of  data  points 

-  Array  of  length  NPT  containing  data  point  abscissae 

-  Array  of  length  NPT  containing  data  point  ordinates 

-  No.  of  knots 

-  Array  of  length  NKT  containing  the  knots 

-  No.  of  integral  ueights  (  -  no.  of  B-splines  - 

order  of  spline  +1  ) 


-  Array  of  length  NUTI  containing  the  integral  ueights 


REAL  X (NPT) .  Y (NPT) , T (NKT) . UTI (NUTI ) , 

*  HL . HR , TL . TR , D1 YL , D1 YR , D2 YDL , D2YDR . D2YTL . D2YTR , 

*  SLOPE  ,Uf1IN,UHAX,UTIAV,  DUMMY 
INTEGER  NPT.NUTI, NKT, K, ID, IT, IU 

K« (NKT -NUT I +1 ) /2 


C  1st  and  2nd  derivatives  at  the  first  tuo  data  points  and  at  the 
C  end-points  of  the  first  knot  interval  are  approximated  by  divided 
C  differences. 


IT-K 

TL-T(K) 

ID-1 

IU-1 

UTI (l)-0.0 
HL-X  (2)  -X  (1) 

HR-X (3) -X (2) 

D1YL-(Y(2)-Y(1))/HL 

01YR- (Y (3) -Y (2) ) /HR 

D2YDL-  (01 YR-D1  YL)  /  (X  (3)  -X  (1 ) ) 

D2YDR-02YDL 

D2YTL-D2YDL 

SLOPE-0.0 
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C  The  next  interval  of  interest  is  the  interval  from  the  right  end  of 
C  the  current  interval  to  the  next  data  point  or  the  next  knot, 

C  uhichever  occurs  first.  The  contribution  to  UTI  from  this  interval 
C  is  determined. 

10  TR-AMIN1 (T (IT+1) ,X (ID+1) ) 

D2YTR«D2YDL+SL0PE*(TR-X (ID) ) 

UTI (IU)-UT1 ( I U) +HL# (D2YTL# (D2YTL+D2YTR) +D2YTR#*2) /3. 

IF (1T.EQ.NUTI+K-11G0  TO  50 
TL-TR 

D2YTI-D2YTR 

IF  (TR.NE.XUD+1) ) GO  TO  30 

ID- ID+1 

D2YDL-D2YDR 

D1YL-D1YR 

HL  -HR 

IF (ID.GT.NPT-2JG0  TO  20 
HR-X ( I D+2) -X ( I D+l ) 

D1 YR-  (Y  (ID+2)  -Y  (ID+1) )  /HR 
D2YDR- (D1YR-01YL) / (HR+HL) 

20  SLOPE- (D2YDR-D2YDD/HL 

30  IF (TR.NE.T (IT+1) ) GO  TO  10 

IT-IT+1 
IU-IU+1 
UTI  (IU)  -0. 0 

40  IF(T(IT+1) .NE.T(IT) )G0  TO  10 

IT-IT+1 
IU-IU+1 
UTI f IU) -0. 0 
GO  TO  40 


C  The  UTI  are  normalized  so  that  most  of  them  are  of  order  1. 

50  UTI AV-SMODAV (NUT1 , UT I ) 

IF(UTIAV.EQ.0.0)GO  TO  70 
UHIN-UTI AV*1.0E-03 
UnAX-UTIAV*1.0E+03 
DO  50  IU-l.NUTI 
DUnriY-UTl (IU) 

IF  ( (UTI  (IU)  .GT.UfllN)  .AND.  (UTI  (IU) .LT.UHAX) )UTI  (IU)  - 
«  UTIAV/UTI  (IU) 

IF (DUnriY.GT.UHAX) UTI  (IUJ-1.0E-03 
50  IF(DUnnv.LE.UniN)UTI (IU)-1.0E+03 

RETURN 

70  DO  80  IU-1.NUTI 

80  UTI  (IU) -1 . 0 

RETURN 
END 
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B.10  XSQC 

FUNCT I  ON  XSQC (N,K, BCOEF , R , VCT , UK . YSQ) 


XSQC  calculates  the  chi-squares 
XSQ-  SITU  (Y(I)-  sum  (BCOEF (J)-Y1(J))*BJ(X(I))  ))**2  *E(1)  ) 

By  subtracting  Y1  from  BCOEF  one  keeps  the  numbers  fairly  small 
thus  avoiding  round-off  error. 

AUTHOR:  David  Hally  ,  May.  1981 

CALLED  by  BSMTH 


REAL  BCOEF (N) . R (K . N) , VCT (N) . UK (N) , XSQC , YSQ 
INTEGER  I.N.K.J.IJl 

XSQC- YSQ 
DO  10  I-l.N 

10  XSQC-XSQC-2. *BC0EF (I )*VCT (I ) 

DO  20  I-l.N 

20  UK  (I) -0.0 

DO  30  I-l.K 

DO  30  J-l.N-I+l 
IJ1-I+J-1 

UK ( J ) -UK ( J) +R ( I , J) *BC0EF (I Jl) 

IFd.EQ.DGO  TO  30 

UK  (I  Jl) -UK  (I  Jl)+R  (I,  J)*BCOEF(J) 

30  CONTINUE 

DO  40  I-1,N 

40  XSQC-XSQC+UK ( I ) *BCOEF ( I ) 

RETURN 

END 
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FIG.  2  SPLINE  CALCULATED  BY  BSMTH:  k  -  4,  STIFFNESS  WEIGHTS  FROM  WTIBEG 
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SPLINE  CALCULATED  BY  BSMTH:  k  -  6,  STIFFNESS  WEIGHTS  FROM  WTIBEG 


FIG.  6  SPLINE  CALCULATED  BY  BSMTH, •  k  -  3  DISCONTINUOUS  DERIVATIVE 
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FIG.  7  SPLINE  CALCULATED  BY  BSMTH  :  k  -  3  DISCONTINUOUS  SPLINE 
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FIG.  9  SOUND  SPEED  PROFILE:  SPLINE  CALCULATED  BY  CUBSPL 
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FIG.  10  SOUND  SPEED  PROFILE:  SPLINE  CALCULATED  BY  3SMTH.  SMOOTHING 
PARAMETER  FROM  PRERR. 
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